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Abstract 

The  work  presented  in  this  thesis  represents  the  latest  in  an  on-going  effort  within 
the  propeller  hydrodynamics  community  to  obtain  rational  design  stage  estimation 
of  dynamically  varying  blade  loads  and  shaft/bearing  forces  for  cavitating  propellers. 
To  this  end,  the  unsteady  cavitating  flow  about  a  marine  propeller  is  treated  in 
nonlinear  theory  by  employing  a  low  order  potential  based  boundary  element  method. 
The  steady  state  oscillatory  solution  is  obtained  by  incremental  stepping  in  the  time 
domain.  The  three  dimensional,  nonlinear  and  time  dependent  boundary  conditions 
are  satisfied  on  an  approximate  boundary  consisting  of  the  propeller  surface  beneath 
the  cavity  and  the  portion  of  the  blade  wake  surface  overlapped  by  the  cavity.  This 
solution  represents  the  first  iteration  of  a  completely  nonlinear  solution  in  which 
the  exact  boundary  conditions  are  satisfied  on  the  exact  flow  boundary,  including 
the  cavity  free-surface.  It  is  then  shown  that  the  convergence  of  such  an  iterative 
solution  for  two-dimensional  flows  is  remarkably  fast  and  that  the  accuracy  of  the 
first  iteration  solution  is  sufficient  for  a  wide  range  of  operating  conditions.  Emphasis 
is  placed  on  developing  an  efficient  and  robust  iterative  scheme  to  predict  general 
cavity  planforms.  The  numerical  method  is  shown  to  be  convergent  and  consistent 
with  fully  nonlinear  results.  Computed  results  are  compared  to  those  from  linear 
theory,  linear  theory  with  leading  edge  corrections,  and  to  published  experimental 
measurements  and  observations. 
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PREFACE 


Cavitation  is  the  process  in  which  vapor  filled  cavities  form  in  a  moving  liquid  in 
response  to  dynamic-pressure  reduction,  as  often  occurs  on  the  back  side  of  propeller 
blades.  Its  importance  to  naval  hydrodynamics  is  paramount,  as  is  perhaps  best 
illustrated  by  the  history  of  the  term  “cavitation.”  The  term  was  coined  (by  Froude) 
during  an  investigation  of  the  failure  of  a  British  destroyer  to  reach  its  design  speed, 
wherein  the  investigators  postulated  that  the  thrust  was  limited  by  cavities  which 
enveloped  the  ship’s  propellers  [41,  62].  The  rest,  as  they  say,  is  history. 

During  the  past  century,  the  study  of  cavitation  has  occupied  both  engineers  and 
mathematicians,  who  tend  to  see  the  world  from  somewhat  different  perspectives. 
Witness,  for  example,  the  views  of  two  men  who  have  made  significant  contributions 
to  the  field:  Professors  R.T.  Knapp  of  California  Institute  of  Technology  and  G. 
Birkhoff  of  Harvard.  Professor  Knapp  (a  civil  engineer),  in  the  preface  to  his  book 
Cavitation  [41],  referred  to  it  as  a  “most  unpleasant”  phenomenon,  with  harmful 
effects  which  “handicap  many  phases  of  science  and  engineering.”  While  he  might 
agree  with  this,  Professor  Birkhoff  (a  mathematician)  looked  more  favorably  on  the 
subject,  as  is  evident  in  the  preface  to  his  book  Jets,  Wakes,  and  Cavities  [4]  where  he 
wrote  of  his  intent  to  report  systematically  on  nearly  a  century  of  “ingenious  research 
in  a  fascinating  and  complex  field.” 

So,  is  cavitation  “unpleasant”,  or  is  it  “fascinating?”  Personally,  I  believe  it 
is  both,  and  I  am  glad  to  have  it  that  way.  On  the  one  hand,  it  is  unpleasant 
enough  to  warrant  a  scientific  investigation,  while  on  the  other  hand  many  of  us  find 
enough  fascination  in  the  topic  to  invest  our  energy,  enthusiasm,  and  wit  in  obtaining 
engineering  solutions.  This  happy  marriage  of  the  practical  and  the  theoretical  is 
true  of  many  engineering  problems  in  naval  hydrodynamics,  a  fact  which  continues 
to  attract  mathematically  inclined  engineers  to  the  field. 
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Nle  number  of  suction  side  panels  on  each  blade  strip  y  plus  the 
number  of  wetted  panels  in  front  of  the  cavity  on  the  pressure 
side 

Nc„  number  of  cavitating  panels  on  the  mth  strip  of  the  wake  surface 
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spanwise  coordinate  normal  to  both  s  and  n  and  forming  a  right- 
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ship  speed 

the  potential  induced  at  the  ith  control  point  on  the  key  blade 
by  a  unit  strength  dipole  at  the  Ith  panel  of  the  mth  strip  of  the 
wake  of  the  kth  blade 

the  potential  induced  at  the  ith  wake  panel  by  a  unit  strength 
dipole  lying  on  the  mth  strip  of  the  wake  of  the  kth  blade 

propeller-fixed  coordinate  system 
ship-fixed  coordinate  system 
hydrofoil  angle  of  attack 

the  openness  of  the  cavity  at  its  trailing  edge,  defined  as  the 
value  of  h  or  hw  at  the  cavity  trailing  edge  (normalized  on  local 
chord  length,  Cm) 

the  tolerance  set  on  8  corresponding  to  a  converged  planform 
perturbation  potential  at  a  field  point  and  an  integration  point, 
respectively 

perturbation  potential  at  the  cavity  leading  edge 


total  potential,  defined  as  the  sum  of  the  inflow  potential  and 
the  perturbation  potential  for  2-D  flows 


$ 


r(r,<)  the  circulation  around  the  blade  at  radius  r  and  time  t 
9  the  circumferential  propeller  coordinate 

9  the  local  angle  between  s  and  v;  s  •  v  =  cos  9 

p  fluid  density 

<r  cavitation  number  based  on  Uoo\  =  f~,?a 

T  OQ 

<rn  cavitation  number  based  on  propeller  rotation;  c  = 

a?  propeller  rotational  speed 


Chapter  1 


Introduction 


A  typical  marine  propeller  operates  in  a  flow  field  which  is  circumferentially  nonuni¬ 
form.  The  variation  in  the  inflow  is  due  primarily  to  the  proximity  of  the  ship’s 
boundary  layer,  but  may  also  be  effected  by  shaft  inclination,  ship  maneuvers,  and 
the  presence  of  upstream  appendages.  This  nonuniformity  often  causes  intermit¬ 
tent  blade  cavitation,  which  in  turn  is  responsible  for  a  host  of  undesireable  effects. 
The  most  conspicuous  of  these  are:  excessive  vibratory  shaft  forces,  unsteady  radi¬ 
ated  pressures  (which  occasionally  result  in  structural  failure  of  nearby  hull  plating), 
blade  surface  erosion,  and  inboard  airborne  noise  levels  in  excess  of  regulatory  stan¬ 
dards.  All  of  these  detrimental  effects  are  associated  with  the  growth  and  collapse  of 
vapor-filled  cavities  and  the  attendant  radiated  pressure  field. 

With  recent  increases  in  the  demand  for  heavily  loaded  efficient  propellers,  cavi¬ 
tation  has  become  less  and  less  avoidable.  As  a  result,  it  has  become  the  task  of  the 
hydrodynamicist  to  predict  the  cavitation  characteristics  of  a  particular  geometry  at 
the  design  stage,  with  the  expectation  that  analytical  methods  can  be  used  to  avoid 
excessive  cavitation  within  an  appropriate  range  of  ship  speed.  Computational  tools 
for  the  reliable  prediction  of  propeller  cavitation  are  therefore  in  high  demand. 

Cavitation  on  marine  propellers  and  hydrofoils  is  a  complex  phenomenon.  It  may 
be  characterized  by  any  combination  of  the  three  forms  (as  delineated  by  Tulin  [66]): 
sheet,  travelling  bubble,  or  cloud  cavitation  (see  Figure  1-1).  A  sheet  cavity  appears 
as  a  single  vapor  bubble  with  a  smooth,  glossy  surface  attached  to  the  blade  near 
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Cavttattng  Tip  Vortex 


Figure  1-1:  A  sketch  of  a  typical  cavitating  propeller,  showing  the  various  types  of 
cavitation. 

its  leading  edge  and  extending  downstream  to  a  region  of  collapse.  The  pressure 
inside  and  on  the  sheet  cavity  is  usually  assumed  to  be  constant  and  equal  to  the 
vapor  pressure.  Although  a  sheet  cavity  can  be  transient,  it  is  often  assumed  that 
the  mechanism  of  unsteadiness  is  the  variations  in  the  inflow  or  in  the  hydrostatic 
pressure.  In  reality,  the  sheet  cavity  is  not  a  smooth,  constant  pressure  surface.  Even 
in  a  uniform  flow  field,  the  sheet  cavity  is  inherently  unsteady  due  to  the  occasional 
(sometimes  periodic)  breakoff  of  the  cavity  trailing  edge  and  the  entrainment  of  trav¬ 
elling  bubbles  in  the  cavity  wake.  Transience  is  also  introduced  at  the  cavity  surface, 
which  in  general  is  a  bubbly,  frothy  region  of  continuous  vaporization  and  condensa¬ 
tion.  Wade  and  Acosta  observed  the  sheet  cavity  on  a  plano-convex  hydrofoil  to  be  a 
rough,  boiling  surface  without  definite  structure  [74].  As  noted  by  Tulin,  however,  it 
is  plausible  that  the  exterior  flow  sees  mainly  the  effect  of  a  coherent  cavity  surface 
on  which  the  pressure  is  constant  [66].  This  observation  is  central  to  justifying  the 
application  of  sheet  cavity  theory,  as  is  done  in  this  thesis. 

In  travelling  bubble  cavitation,  distinct  vapor  bubbles  grow  and  collapse  over 
the  surface  of  the  blade  as  they  travel  from  regions  of  low  pressure  to  regions  of  high 
pressure.  The  coalescence  of  clusters  of  travelling  bubbles  seems  to  play  an  important 
role  in  the  formation  of  a  sheet  cavity  in  some  circumstances.  However,  this  process 
is  not  well  understood.  Modelling  bubble  cavitation  is  also  problematic.  While  the 
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growth  and  collapse  of  a  single  spherical  bubble  is  treated  by  Rayleigh- Plesset  theory 
[41],  no  adequate  theory  exists  to  treat  clusters  of  individual  bubbles  in  the  presence 
of  solid  boundaries.  Nonetheless,  some  progress  has  been  made  in  recent  years  in 
applying  analytical  and  numerical  methods  to  model  the  behavior  of  multiple  bubble 
clusters  [6,  8,  19).  An  example  is  the  work  of  Chahine,  who  has  used  a  boundary- value 
approach  and  applied  a  surface  panel  method  to  multiple  collapsing  bubbles  [6], 

Downstream  of  a  collapsing  sheet  cavity  is  often  observed  a  wake-entrained  cluster, 
or  cloud,  of  very  small  bubbles.  The  collapse  of  the  cavitation  cloud  —  usually  near 
the  trailing  edge  of  a  propeller  blade  where  the  cloud  is  subject  to  a  positive  pressure 
gradient —  is  known  to  cause  severe  surface  erosion.  The  phenomenon  has  thus  far 
defied  mathematical  modelling. 

Associated  with  cavitating  propellers  are  hub  and  tip  vortices  which  may  have 
vapor-filled  cores.  The  cavitating  tip  vortex  often  exceeds  in  volume  the  sheet  cavity 
and  has  a  marked  effect  on  the  local  flow  at  the  tip.  The  process  of  cavitation 
inception  is  not  well  understood,  nor  are  the  roles  of  viscosity  or  nucleation.  Propeller 
sheet  cavities  in  nonuniform  flow  are  sometimes  observed  to  collapse  inward  from  the 
leading  edge  and  to  disintegrate  into  a  cloud  cavity.  The  effect  of  viscosity  on  the 
cavity  flow,  which  has  been  observed  mainly  to  be  a  determinant  of  the  point  of 
detachment  [2,  13],  has  yet  to  be  included  in  any  propeller  analysis  tool. 

At  the  present  time,  there  are  essentially  two  distinct  avenues  for  modelling  cav¬ 
itation.  One  is  to  treat  the  cavity  strictly  as  a  constant  pressure  sheet,  while  the 
other  is  to  consider  the  cavity  to  be  a  cluster  of  bubbles.  An  example  of  the  latter 
approach  is  the  work  of  Kubota  et  al  [43],  who  solve  the  compressible  Navier-Stokes 
equations  while  treating  the  two-phase  flow  as  a  single-phase  compressible  fluid  with 
drastically  varying  density.  They  allow  a  bubble  cluster  to  grow  and  collapse  accord¬ 
ing  to  Rayleigh’s  pressure/ volume  relationship  for  noninteracting  spherical  bubbles. 
Initial  results  compare  qualitatively  well  with  experiments.  However,  the  method  is 
hampered  by  a  bubble  model  which  does  not  capture  the  effects  of  bubble-bubble 
interaction  or  those  of  bubble  non-axisymmetry.  It  has  been  shown  that  these  ef¬ 
fects  are  very  important  to  the  prediction  of  bubble  collapse  and  the  resulting  high 
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impact  point  loads  and  generation  of  acoustic  noise  [6].  Nonetheless,  with  inevitable 
improvements,  this  method  may  lead  the  way  to  a  more  complete  understanding  of 
the  phenomenon. 

The  approach  taken  in  this  work  is  to  treat  propeller  cavitation  strictly  as  sheet 
cavitation.  The  advantage  of  this  approach  lies  in  its  relative  mathematical  simplic¬ 
ity.  For  example,  if  the  flow  is  assumed  to  be  potential  in  nature  (irrotational  and 
incompressible  flow  of  an  inviscid  fluid),  then  the  governing  equation  in  steady  flow 
is  a  linear  elliptic  partial  differential  equation  for  a  single  scalar.  The  use  of  surface 
singularities  and  the  now  classic  application  of  Green’s  3rd  identity  are  well  suited 
for  such  a  problem.  The  disadvantages  of  this  approach  are  the  neglect  of  viscosity, 
vortex  cavitation,  bubble  cavitation,  and  cloud  cavitation.  However,  there  are  two 
additional  rationales  for  taking  the  potential  flow  avenue.  First,  it  is  believed  that  the 
first-order  contributor  to  dynamically  varying  blade  loads  and  radiated  pressures  is 
the  transient  sheet  cavity  [66].  Second,  many  or  all  of  the  neglected  phenomena  may 
be  included  as  refinements  to  the  potential  flow  solution.  For  example,  a  viscous- 
inviscid  boundary  layer  interaction  could  be  used  to  model  the  effects  of  viscosity. 
As  another  example,  a  local  tip  solution,  including  a  cavitating  tip  vortex,  could  be 
matched  to  the  outer  sheet  cavity  solution.  The  potential  of  such  a  model  to  cap¬ 
ture  most  of  the  physics  at  reasonable  computational  expense  is  the  motivation  for 
obtaining  an  accurate ,  efficient,  and  robust  potential  flow  solution.  This  is  the  goal 
of  this  thesis. 


1.1  Objectives 

The  objective  is  to  develop  a  robust  and  computationally  efficient  nonlinear  method 
of  predicting  unsteady  propeller  sheet  cavitation  and  its  resulting  free- space  potential 
field.  The  approach  is  to  use  a  low  order  potential  based  boundary  element  method 
to  obtain  the  steady  state  oscillatory  sheet  cavity  solution,  including  an  accurate 
prediction  of  the  cavity  history.  Emphasis  is  placed  on  obtaining  accurate  prediction 
of  cavity  planforms  and  shapes.  It  is  expected  that  such  a  solution  will  improve  the 
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accuracy  of  the  prediction  of  cavitation-induced  unsteady  blade  load  distributions  and 
hull  vibratory  excitations.  This  in  turn  will  improve  the  identification  of  unacceptable 
propeller  or  hull  designs  early  in  the  design  process,  resulting  in  a  more  cost-effective 
design. 


1.2  Fundamental  Assumptions 

The  propeller  blades  are  assumed  to  be  rigid  bodies  arranged  symmetrically  around  a 
common  axis  and  rotating  with  constant  angular  velocity  in  an  unbounded  fluid.  The 
presence  of  other  solid  boundaries,  such  as  the  hull,  duct  or  struts,  is  ignored.  The 
hub,  however,  may  be  included  in  the  analysis.  The  inflow  is  known  and  represents 
the  wake  field  of  the  ship.  It  is  assumed  to  be  the  effective  wake  which  includes  the 
interactions  between  the  vorticity  of  the  inflow  in  the  absence  of  the  propeller  (the 
nominal  wake)  and  the  vorticity  due  to  the  propeller  [29,  27].  The  inflow  is  assumed 
to  be  constant  over  the  axial  extent  of  the  propeller. 

The  fluid  is  assumed  to  be  inviscid  and  free  of  impurities  and  the  resulting  flow  to 
be  incompressible  and  irrotational.  The  cavity  is  assumed  to  be  a  constant-pressure 
surface  which  grows  on  the  suction  side  of  a  blade  when  the  local  pressure  falls  below 
the  vapor  pressure  and  collapses  when  it  rises  above  it.  The  pressure  side  of  the  blade 
is  assumed  to  be  cavity-free.  The  cavity  is  assumed  to  be  a  sheet  cavity;  neither 
travelling  bubble  nor  cloud  cavitation  is  included  in  the  model.  Cavitating  tip  and 
hub  vortices  are  also  excluded  from  the  analysis.  The  detachment  line  of  the  sheet 
cavity  (the  locus  of  all  spanwise  points  where  the  cavity  begins)  is  assumed  to  be 
known. 

Despite  the  drastic  simplification,  the  problem  at  hand  remains  a  daunting  one. 
The  main  difficulty  in  the  analysis  of  sheet  cavitation  arises  from  the  need  to  de¬ 
termine  the  free  streamline  (the  cavity  surface)  on  which  the  pressure  is  prescribed. 
The  fact  that  the  location  of  a  portion  of  the  flow  boundary  is  unknown  makes  the 
problem  strongly  nonlinear.  The  inherent  complexity  of  three  dimensional  flows  adds 
to  the  difficulty  in  the  analysis.  Variation  in  cavity  extent  is  often  nonlinear  in  time, 
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since  a  typical  propeller  blade  operating  in  a  ship  wake  goes  from  fully  wetted  to 
supercavitating  and  back  to  fully  wetted  within  a  fraction  of  a  revolution. 

Perhaps  it  was  these  compounded  complexities  which  prompted  Batchelor  to  state 
the  following: 

If  the  shape  of  the  permanent  cavity  under  different  conditions  could  be 
predicted,  the  propeller  could  be  designed  to  allow  for  the  existence  of  the 
cavity,  but  this  is  not  possible  in  general. 

G.K.  Batchelor,  in  An  Introduction  to  Fluid  Dynamics  [3] 

Batchelor  was  only  partly  correct  in  this  statement.  Propeller  design  is  still  done  with 
only  approximate  account  taken  of  the  effects  of  cavitation.  However,  advances  in 
high  speed  computing  have  allowed  for  the  analysis  of  three  dimensional  sheet  cavities 
on  given  geometries.  Prediction  of  propeller  sheet  cavitation  may  now  be  included  in 
the  overall  design  process. 


1.3  Previous  Research 

1.3.1  2-D  Free  Streamline  Flows 

The  solution  of  the  2-D  free  streamline  problem  began  with  Helmholtz  and  Kirchoff, 
who  solved  for  the  flow  around  flat  plate  or  polygonal  bodies  at  zero  cavitation  num¬ 
ber.  These  were  exact  solutions  which  made  use  of  the  hodograph  technique,  which 
was  introduced  by  Helmholtz  more  than  a  century  ago  [4].  The  method  was  later 
extended  to  treat  curved  bodies  at  zero  cavitation  number  by  Levi-Civita  [4].  The 
theory  was  applied  by  Wu  and  Wang  [76]  and  later  by  Furuya  [14]  for  the  analysis  of 
supercavitating  hydrofoils  in  the  presence  of  a  free  surface. 

The  extension  to  non-zero  cavitation  number  resulted  in  a  diversity  of  termination 
models.  Among  these,  the  most  popular  are  the  Riabouchinsky  model  [58, 60,  61],  the 
re-entrant  jet  model  [9,  42],  the  open  wake  model  [10],  and  the  spiral  vortex  models 
[65].  Among  researchers,  there  has  been  no  consensus  regarding  which  model  is  best 
suited  to  represent  cavity  flows.  A  comprehensive  review  of  termination  models  may 
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be  found  in  [75]  and  [63].  A  description  of  termination  models  developed  in  this  thesis 
is  included  in  Appendix  C. 

Despite  producing  exact  solutions,  the  hodograph  method  with  the  various  ter¬ 
mination  models  were  never  widely  used  due  to  the  inherent  difficulty  in  treating 
general  shaped  bodies  and  the  restriction  to  planar  flows.  These  deficiencies  led  to 
the  introduction  of  linear  theory. 

1.3.2  2-D  Linear  Theory 

The  solution  of  two-dimensional  cavity  flows  for  general  geometries  became  more  ac¬ 
cessible  with  the  development  of  linear  theory,  as  introduced  by  Tulin  for  general 
shape  supercavitating  hydrofoils  at  zero  cavitation  number  [63],  and  later  at  non-zero 
cavitation  number  [64].  Linear  theory  was  first  applied  to  partially  cavitating  hydro¬ 
foils  by  Acosta  [1]  and  Geurst  and  Timman  [16],  independently.  The  linear  cavity 
theory,  which  closely  resembles  the  classical  thin- wing  theory,  assumes  thin  cavity  and 
foil  thickness  relative  to  the  foil  chord  length  and,  as  a  result,  simplifies  the  dynamic 
boundary  condition  from  the  requirement  of  constant  total  velocity  to  that  of  con¬ 
stant  horizontal  perturbation  velocity.  This  boundary  condition  is  applied,  together 
with  the  linearized  wetted  surface  kinematic  boundary  condition,  on  a  projection  of 
the  hydrofoil  surface  on  the  free  stream  axis.  Linear  cavity  theory  has  been  applied 
by  several  researchers  for  the  analysis  of  partial  and  supercavitating  flows  around 
special  and  general  shape  2-D  hydrofoil  geometries.  A  list  of  references  which  address 
the  application  of  linear  cavity  theory  may  be  found  in  [67]  or  [33]. 

1.3.3  3-D  Cavity  Flows 

The  flow  around  cavitating  finite  span  hydrofoils1  was  first  treated  in  a  stripwise  sense 
with  the  three  dimensional  flow  effects  introduced  by  matching  the  inner  solution 
(from  either  linear  or  nonlinear  theory)  with  the  solution  from  lifting  line  theory  in 
the  outer  domain.  This  approach  was  used  to  study  supercavitating  hydrofoils  by 

tAIso  called  3-D  hydrofoils  in  this  work. 
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Nishiyama  [56],  Leehey  [48]  and  Furuya  [15]  and  partially  cavitating  hydrofoils  by 
Uhlman  [68].  However,  this  method  is  valid  only  for  high  aspect  ratio  hydrofoils  and 
its  accuracy  deteriorates  near  the  hydrofoil  tips.  The  complete  three  dimensioned 
flow  effects,  in  the  context  of  linear  cavity  theory,  were  included  in  the  analysis  of 
supercavitating  3-D  hydrofoils  by  Widnall  [59],  who  employed  a  numerical  pressure 
doublet  and  source  lifting  surface  technique.  Jiang  and  Leehey  instead  applied  a 
discrete  vortex  and  source  lattice  lifting  surface  technique  [25] ,  and  also  introduced  an 
iterative  scheme  to  determine  a  cavity  planform  corresponding  to  a  uniform  cavitation 
number  spanwise.  The  use  of  the  vortex/source  lattice  led  to  the  treatment  of  fully 
wetted  propellers  in  unsteady  flow  by  Kerwin  and  Lee  [32],  which  was  then  extended 
to  the  analysis  of  unsteady  sheet  cavitation  by  C-S.  Lee  [45]  and  Breslin  et  al  [5]. 
These  researchers  also  introduced  the  ability  to  treat  mixed  cavity  patterns  in  partial 
and  supercavitation. 

1.3.4  Nonlinear  Solutions 

Due  to  the  success  of  thin  wing  theory  in  the  analysis  of  fully  wetted  flows,  one  might 
assume  that  linear  cavity  theory  would  be  appropriate  for  analyzing  hydrofoil  or  pro¬ 
peller  sheet  cavitation.  However,  a  serious  defect  of  the  linear  theory  arises  in  the 
case  of  partially  cavitating  flows  around  hydrofoils  with  rounded  leading  edges.  In 
particular,  linear  theory  predicts  that  the  cavity  extent  and  volume  should  increase 
with  increasing  thickness  provided  that  the  flow  conditions  remain  the  same.  On  the 
contrary,  the  short  cavity  theory  of  Tulin  and  Hsu  [67]  and  the  numerical  nonlin¬ 
ear  surface  vorticity  method  developed  by  Uhlman  [70],  predict  that  the  cavity  size 
should  decrease  with  increasing  thickness.  In  addition,  it  is  known  from  experimental 
evidence  that  increasing  the  leading  edge  radius  delays  or  reduces  leading  edge  cavi¬ 
tation.  This  defect  in  the  linear  cavity  theory  was  recently  addressed  by  Kinnas  [33], 
who  introduced  the  nonlinear  leading  edge  correction  to  account  for  the  breakdown  of 
the  linear  cavity  solution  in  the  vicinity  of  a  round  leading  edge.  The  nonlinear  blade 
thickness  effects  have  also  been  included  in  the  linearized  method  for  the  analysis  of 
unsteady  propeller  cavitation  [45]  by  implementing  the  leading  edge  correction  locally 
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at  planes  normal  to  the  blade  leading  edge  [31,  34]. 

Numerical  nonlinear  solutions  have  evolved  because  of  the  defect  in  the  linear  the¬ 
ory.  Various  velocity-based  boundary  element  methods  (BEM’s)  have  been  applied 
for  the  analysis  of  the  flow  around  2-D  and  3-D  partially  and  supercavitating  hydro¬ 
foils.  These  include  methods  by  Pellone  and  Rowe  [57],  for  supercavitating  2-D  and 
3-D  hydrofoils;  by  Uhlman  [70,  71],  for  partially  and  supercavitating  2-D  hydrofoils; 
and  by  Lemonnier  and  Rowe  [20],  for  partially  cavitating  2-D  hydrofoils.  In  each 
of  these  BEM’s,  the  exact  cavity  boundary  is  determined  by  an  iterative  process  in 
which  the  dynamic  boundary  condition  is  satisfied  on  an  approximate  cavity  surface 
and  the  kinematic  boundary  condition  is  used  to  update  the  surface.  The  iterations 
terminate  when  both  the  kinematic  and  the  dynamic  boundary  conditions  are  sat¬ 
isfied  on  the  computed  cavity  surface.  In  order  to  solve  the  problem,  all  of  these 
techniques  assume  that  the  cavity  extent  is  known;  the  cavitation  number  is  part  of 
the  solution. 

Each  of  the  numerical  nonlinear  methods  requires  modest  computer  effort  to  solve 
the  2-D  problem,  but  their  extension  to  3-D  was  limited  due  to  slow  rates  of  con¬ 
vergence  to  the  fully  nonlinear  solution.  If  a  numerical  nonlinear  method  was  to  be 
applied  to  3-D  geometries,  a  better  tool  was  needed. 

1.3.5  Present  Method 

In  the  present  method,  a  potential  based  boundary  element  method  (t.e.,  one  based 
on  Green’s  third  identity  for  the  perturbation  potential)  is  applied  to  the  propeller 
sheet  cavity  problem  in  unsteady  flow.  The  method  was  first  applied  to  find  the  fully 
nonlinear  fixed-length  solution2  for  partially  and  supercavitating  2-D  hydrofoils  in 
steady  flow,  as  reported  in  [35].  An  iterative  method  was  employed  to  find  the  exact 
cavity  shape,  similar  to  the  one  described  in  the  previous  section  for  the  velocity- 
based  methods.  In  the  first  iteration,  the  panels  representing  the  cavity  were  placed 
on  the  foil  surface.  In  subsequent  iterations,  the  cavity  surface  was  updated  and  the 

JIn  the  fixed-length  solution,  the  length  is  specified  and  the  cavitation  number  is  unknown. 
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First  iteration  panelling 


Subsequent  iteration  panelling 


Figure  1-2:  Panel  arrangements  for  the  first  iteration  (left)  and  subsequent  iterations 
(right). 

panels  moved  to  the  updated  sui.'ace  (see  Figure  1-2).  At  that  time,  the  method 
was  characterized  by  very  rapid  convergence  to  the  fully  nonlinear  solution  and  the 
remarkable  accuracy  of  the  first  iteration.  The  convergence  of  the  cavitation  number 
and  cavity  volume  with  iterations  is  much  faster  for  the  potential  based  method  than 
for  the  velocity  based  methods.  This  contrast  is  shown  in  Figure  1-3,  which  shows 
the  results  of  the  present  potential  based  method  and  the  results  of  Uhlman’s  velocity 
based  method  [71]  for  a  NACA16006  hydrofoil  at  angle  of  attack  a  =  4°  and  cavity 
length-to-chord  ratio  ~  =  0.5. 

In  a  similar  method  developed  at  the  same  time  by  Lee  et  al  [47],  similar  behavior 
of  the  solution  with  iterations  was  found.  However,  in  that  work  a  termination  model 
was  not  applied.  The  importance  of  applying  a  termination  model  in  the  nonlinear 
solution  is  discussed  in  Appendix  C. 

In  addition  to  the  quick  convergence,  the  first  iteration  cavity  area  differed  from 
the  converged  nonlinear  area  by  no  more  than  2%  (see  Appendix  B)  for  partially 
cavitating  hydrofoils  at  a  wide  range  of  angles  of  attack.  It  was  recognized  then  that 
such  accuracy  could  represent  substantial  savings  and  render  a  fully  three  dimensional 
nonlinear  solution  computationally  feasible  if  a  single  iteration  toward  the  nonlinear 
solution  is  sufficient. 

The  method  was  then  extended  to  treat  2-D  and  3-D  partially  cavitating  hydrofoils 
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Figure  1-3:  Convergence  with  number  of  iterations  of  the  cavitation  number  and  the 
cavity  volume  for  a  NACA16006  foil  at  a  =  4°,  I  =  .5.  Results  are  shown  for  the 
present  method  (left)  and  Uhlman’s  method  (right). 

for  the  more  reuli  .ic  case  where  the  cavitation  number  is  known  and  the  cavity 
extent  must  oe  determined  [37].  In  the  3-D  extension,  the  fully  three  dimensional, 
nonlinear  and  time  dependent  boundary  conditions  were  applied.  The  development  of 
an  efficient  and  robust  iterative  algorithm  to  find  the  cavity  extent  was  central  to  this 
work.  Next,  the  model  was  extended  to  treat  supercavitation,  along  with  the  more 
difficult  problem  of  treating  mixed  partial/ supercavity  patterns.  The  emphasis  at  this 
stage  was  to  develop  a  uniform  analytical  approach,  one  that  made  use  of  a  single 
algorithm  to  compute  general  shape  planforms  including  partial  and  supercavitation. 
The  method  was  tested  extensi  /ely  for  convergence  and  consistency.  The  results  of 
this  part  of  the  research  were  reported  in  [13]. 

In  the  3-D  hydrofoil  nonlinear  solution,  only  the  first  iteration  is  performed.  Based 
on  the  accuracy  of  the  2-D  first  iteration  solution,  it  was  deemed  unnecessary  to  con¬ 
tinue  iterating.  Obviously,  this  decision  was  also  motivated  by  the  extensive  computer 
effort  necessary  in  solving  3-D  flows.  Despite  the  approximation,  the  nonlinear  char¬ 
acter  of  the  solution  has  been  preserved,  and  the  solution  is  found  with  moderate 
CPU  expense.  In  the  remainder  of  this  thesis,  the  solution  will  be  called  nonlinear 
in  order  to  differentiate  it  from  the  linear  solution.  The  converged  nonlinear  solution 
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will  be  referred  to  as  fully  nonlinear. 

The  latest  extension  of  the  method  involved  the  treatment  of  cavitating  propellers, 
in  either  uniform  or  nonuniform  inflow.  The  details  of  the  method  will  be  described 
in  this  thesis  mainly  in  the  context  of  propellers  in  nonuniform  flow. 
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Chapter  2 


Mathematical  Formulation 


In  this  chapter,  a  precise  definition  of  the  problem  is  written  down.  Green’s  formula 
is  applied  on  an  approximate  flow  boundary,  which  includes  a  portion  of  the  zero¬ 
thickness  wake  sheet,  and  this  necessitates  special  treatment  of  the  boundary  integral 
equation  in  order  to  extract  the  correct  local  contributions.  Following  this,  detailed 
derivations  of  the  fully  three-dimensional  nonlinear  boundary  conditions  are  provided. 

2.1  Definition  of  the  Problem 

Consider  a  cavitating  propeller  subject  to  a  nonuniform  inflow  Uw{zs,rs,9s),  with 
the  subscript  5  denoting  the  ship-fixed  coordinate  system  in  which  the  wake  is  defined. 
As  mentioned,  Uw  is  assumed  to  be  the  effective  wake.  Some  details  of  the  geometry 
are  shown  in  Figure  2-1.  The  solution  is  found  in  the  (x,  y,  z)  coordinate  system,  which 
rotates  with  the  propeller.  The  propeller  is  assumed  (without  loss  of  generality)  to 
be  right-handed  and  to  rotate  at  a  constant  angular  velocity  u).  The  inflow  relative 
to  the  propeller  is 


Uin(x,y,z,t)  =  Uw{x,r,9  -ut)  +  u>  x  *  (2.1) 

where  r  =  s/y2  +  z2,  9  =  arctan(^),  and  a  =  (x,y,z).  The  inflow  velocity  Uw  is 
expressed  in  terms  of  its  axial,  radial  and  tangential  components  (labeled  Vm,  Vr,  and 
Vt,  respectively).  The  components  are  defined  in  terms  of  harmonic  coefficients  as 
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follows 


KM) 

Vs 

KM) 

Vs 

KM) 

Vs 


NHARMA  NHARMA 

A*(r)  +  4*(r)  cos  ^  B*(r)  sin  nB 

n=2  n=2 

NHARMR  NHARMR 

j4i(r)+  5Z  -^n(r)  cos  nB  +  ^  B„(r)sinnB 

n= 2  n=2 

NHARMT  NHARMT 

A‘i(r)  4-  -^n(r)  cos  nB  +  ^  B^(r)sinn0  (2.2) 

n=2  n=2 


where  .An  and  J3„  are  the  harmonic  coefficients  and  the  superscripts  x,r  and  t  denote 
axial,  radial  and  tangential  components  [28].  The  leading  coefficient  A\  is  the  mean 
value.  The  angular  coordinate  B  is  positive  in  the  clockwise  direction  looking  down¬ 
stream  (see  Figure  2-1).  NHARMA ,  NHARMR  and  NHARMT  are  the  number  of 
harmonics  which  are  used  to  describe  a  given  inflow.  The  inflow  velocity  components 
vary  only  radially  and  circumferentially;  they  are  assumed  to  be  constant  over  the 
axial  extent  of  the  propeller. 

For  the  moment,  assume  that  the  propeller  has  a  developed  sheet  cavity  whose 
time-dependent  surface  is  denoted  by  Sc  (*),  as  shown  in  Figure  2-1.  The  fluid  is 
assumed  to  be  inviscid  and  the  resulting  flow  to  be  incompressible  and  irrotational. 
In  this  case,  the  time-dependent  total  flow  velocity  q(x,y,z,t ),  can  be  written  in 
terms  of  the  perturbation  potential,  (j>{x,y,z,t),  as  follows: 


q{x,y,z,t)  =  Uin(x,y,z,t)  +  V<^(x,y,z,t).  (2.3) 

The  goal  is  to  determine  the  potential  field  <f>(x,y,z,t),  as  well  as  the  cavity  surface 
Sc{t).  Once  <f>  is  known,  the  pressure  distribution  may  be  computed  by  numerically 
differentiating  the  potentials  and  applying  Bernoulli’s  equation.  The  unsteady  blade 
load  distribution  may  then  be  determined  by  integrating  the  pressure.  Knowing  Sc  (0> 
the  cavity  volume  history  may  be  computed.  In  the  following  sections,  the  equations 
and  boundary  conditions  necessary  for  the  solution  will  be  derived. 
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Figure  2-1:  The  propeller  and  cavity  geometry,  coordinate  systems,  and  nonuniform 
inflow.  xs,ys,zs  represent  the  ship-fixed  coordinate  system;  x,y,  z  represent  the 
propeller-fixed  coordinate  system.  The  point  p  is  a  point  on  the  propeller  surface, 
defined  by  the  vector  x  from  the  origin  of  the  propeller-fixed  system. 
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Figure  2-2:  Definition  of  the  flow  boundaries  on  which  the  boundary  conditions  should 
be  satisfied. 

2.2  Green’s  Formula 

The  perturbation  potential  <f>p(x,  y,  z,  t)  at  any  point  p  which  lies  either  on  the  wetted 
blade  (or  hub)  surface1,  Sws(t )>  or  on  the  cavity  surface,  Sc{t),  (both  surfaces  are 
shown  in  Figure  2-2)  must  satisfy  Green’s  third  identity: 


2*  4>p(t) 


=  / 

Snrs(t)uSc(t) 


9G(r,q) 

dnq(t) 


G(p;9)|^U) 


dS  + 


+ 


J  A<j>w(rq}9q,t) 

Sw(t) 


dnq(t) 


dS\ 


p  €  ( Sws  U  Sc) 


(2.4) 


The  subscript  q  corresponds  to  the  variable  point  in  the  integrations.  nq(t)  is  the  unit 
vector  normal  to  the  surface  of  the  propeller,  the  cavity,  or  the  wake.  The  unit  normal 
points  into  the  fluid  (on  the  wake  surface,  nq(t)  is  oriented  such  that  it  points  in  the 
same  general  direction  as  the  normal  on  the  auction  side  of  the  blade).  A ^„(r,0,f)  is 
the  potential  jump  across  the  wake  sheet,  Sw(t),  and  G(p;q)  =  1  /R(p',q)  is  Green’s 
function,  where  R(p;  q)  is  the  distance  from  the  field  point  p  to  the  variable  point  q. 

Equation  (2.4)  expresses  the  perturbation  potential  on  the  surface  formed  by  the 

1The  wetted  blade  surface  is  that  part  of  the  blade  which  is  not  eavit sting. 
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union  of  the  cavity  and  blade  surfaces,  Sws(t)  U  Sc(t),  as  a  superposition  of  the 
potentials  induced  by  a  continuous  source  distribution,  G ,  and  a  continuous  dipole 
distribution,  on  Sws(t)  U  Sc(0>  and  a  continuous  dipole  distribution  on  the 
trailing  wake  surface,  Sw(t).  The  application  of  Green’s  3r<^  identity  to  problems  in 
potential  flow  is  classic  [44,  50,  54],  while  applications  to  propeller  and  rotor  flows 
are  more  recent  [49,  51,  21,  22,  30,  38]. 

Ultimately,  the  exact  nonlinear  potential  flow  solution  will  be  found  when  the 
kinematic  and  dynamic  boundary  conditions  (which  may  be  applied  simultaneously 
with  Green’s  formula  (2.4))  are  satisfied  on  the  exact  flow  boundary.  However,  we 
face  the  usual  problem  that  the  position  of  the  cavity  surface  is  unknown.  As  a  first 
iteration  towards  the  fully  nonlinear  solution,  we  may  apply  the  boundary  conditions 
on  an  approximate  flow  boundary.  A  natural  choice  of  an  approximate  boundary  is 
defined  as  follows: 

•  it  coincides  with  the  blade  surface  (including  the  part  of  the  blade  beneath  the 
cavity).  This  surface  will  be  referred  to  as  Scb • 

•  if  the  blade  is  supercavitating,  the  upper  and  lower  sides  of  the  supercavity 
downstream  of  the  blade  trailing  edge  are  collapsed  to  a  single  surface.  The 
two  sides  of  the  collapsed  cavity  surface  coincide  with  the  two  sides  of  the  zero- 
thickness  trailing  wake  sheet  (see  discussion  below).  This  surface  will  be  referred 
to  as  5cw(0- 

A  sketch  of  the  approximate  boundary  is  shown  in  Figure  2-3.  The  geometry  of 
the  wake,  Sw(t)  is  assumed  to  be  invariant  in  time  and  is  taken  to  be  the  same  as 
the  steady-flow  relaxed  wake  corresponding  to  the  circumferentially  averaged  inflow 
[18].  The  approximate  flow  boundary  therefore  coincides  with  that  of  the  fully  wetted 
solution,  as  described  by  Hsin  [22]. 

As  a  result  of  the  linearization  of  the  part  of  the  supercavity  downstream  of  the 
blade  trailing  edge,  the  surface  Scw(t)  may  be  considered  to  be  a  force-free  surface. 
The  pressure  across  the  collapsed  cavity  surface  must  therefore  be  continuous  and 
equal  to  the  cavity  pressure  pe.  The  pressure  across  the  blade  wake  surface  must  also 
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Figure  2-3:  The  approximate  cavity  surface  on  which  the  boundary  conditions  are 
applied. 

be  continuous.  Therefore,  it  is  possible  to  consider  these  two  surfaces  to  coincide, 
with  the  condition  of  pressure  continuity  being 

P+  =  P~  =  Pc  on  Sew  (t) 

P+  =  p~  =p  on  Sw(t).  (2.5) 

As  mentioned  in  Chapter  1,  satisfying  the  boundary  conditions  on  the  approxi¬ 
mate  boundary  serves  as  the  first  iteration  in  a  fully  nonlinear  solution.  In  the  fully 
nonlinear  solution,  subsequent  iterations  are  found  by  satisfying  the  dynamic  bound¬ 
ary  condition  on  an  updated  cavity  surface  (the  kinematic  boundary  condition  is  used 
to  update  the  surface,  as  described  in  section  2.4).  The  solution  is  then  considered 
converged  when  the  cavity  surface  does  not  change  (to  within  a  tolerance)  between 
two  consecutive  iterations.  It  was  also  mentioned  in  Chapter  1  that,  using  the  present 
potential  based  panel  method,  the  first  iteration  solution  (where  Green’s  formula  is 
satisfied  on  the  approximate  boundary)  is  extremely  close  to  the  converged  nonlinear 
solution  for  a  wide  range  of  operating  conditions.  As  a  result,  it  is  deemed  unnec¬ 
essary  to  go  beyond  the  first  iteration  towards  the  fully  nonlinear  solution.  In  view 
of  the  high  computational  cost  involved  in  regridding  and  re-computing  influence 
coefficients,  the  importance  of  this  conclusion  cannot  be  overstated. 
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There  is  a  simple  recurring  theme  throughout  this  thesis:  avoid  recomputing  in¬ 
fluence  coefficients.  In  fact,  the  struggle  to  avoid  computing  additional  influence 
coefficients  is  the  driving  force  behind  many  of  the  numerical  techniques  to  be  in¬ 
troduced  in  this  thesis.  In  the  fully  wetted  steady  flow  solution,  the  computation  of 
influence  coefficients  dominates  the  computation  time.  The  solution  of  the  unsteady 
cavitating  propeller  will  be  found  in  the  time  domain,  with  the  fully  wetted  solution 
serving  as  the  initial  condition.  A  significant  amount  of  computation  is  avoided  by 
using  the  same  influence  coefficients  in  both  the  fully  wetted  and  cavitating  solu¬ 
tions.  A  computationally  efficient  solution  is  made  feasible  by  the  accuracy  of  the 
first  iteration. 

Green’s  formula  (2.4)  holds  on  the  exact  flow  boundary.  It  does  not,  however, 
apply  on  the  upper  surface  of  the  wake,  which  is  part  of  the  approximate  boundary. 
Care  must  be  taken  in  extracting  the  correct  local  contribution  when  the  field  point 
lies  on  the  wake  surface. 

Considering  the  approximate  boundary,  the  first  integral  on  the  right-hand-side 
of  (2.4)  may  be  decomposed  into  integrals  over  the  blade  surface  and  the  portion 
of  the  wake  which  is  overlapped  by  the  cavity  (previously  denoted  Scb  and  Sew , 
respectively,  in  Figure  2-3).  The  exact  form  of  Green’s  formula  depends  on  the 
location  of  the  field  point,  which  will  either  be  on  Scb  or  on  Sew •  Each  case  will  be 
considered  separately. 

Field  Point  on  Scb 

If  the  field  point  is  on  the  blade  surface  Scb,  the  local  contribution  is  extracted  from 
the  first  integral  in  (2.4)  and  Green’s  formula  becomes 

-  /  °M£+(,>-5W"  + 


41 


(2.6) 


+  j  (A <f>w(r,0,t))dGQ^-dS 

ScwuSh' 

where  the  superscripts  -f  and  -  correspond  to  the  upper  and  lower  sides  of  the  wake 
surface,  respectively.  Note  that  the  wake  surface  has  been  divided  into  Sew  and  Sw, 
as  shown  in  Figure  2-3. 

The  velocity  normal  to  Sew  is  discontinuous  across  the  surface.  The  jump  in 
defines  a  source  distribution,  of  density  qw{t),  which  represents  the  cavity  thickness: 


ftyj  i*\  _  (t) 

-  Qn  (<) 


dna 


(2.7) 


The  potential  is  also  discontinuous  across  the  Sw,  where  the  jump  A <f>w(r,0,t)  is 
related  to  the  local  circulation  history.  This  is  described  in  section  4.4. 

Inserting  (2.7)  into  (2.6)  yields: 


ScB 


9G{p\  g) 

0nq 


dS~ 


J  qw{t)G{p;  q)  dS  +  J  A<i>w(r,0,t) — ^  —dS  on  SCb- 

Scw(t)  ScwuSw 

(2.8) 


Field  Point  on  Sew 

Consider  the  case  where  the  field  point  is  on  the  upper  side  of  the  collapsed  cavity 
surface  in  the  wake,  SCw(t).  In  this  case,  local  contribution  must  be  extracted  from 
the  integral  over  Sw(t)  on  the  right  hand  side  of  (2.4),  which  reduces  to 


SCB  L  ’  ’ 


dS  - 


~  J  <lw{t)G(p’,q)  dS  +  j  A<j>w(r,0,t) — -d5  on  SCw- 

Scw(t)  ScwuSw 


(2.9) 
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Derivation  of  this  result  is  provided  in  Appendix  A. 

The  left  hand  side  of  (2.9)  may  be  written  in  terms  of  the  potential  on  one  side 
of  the  wake  sheet  and  the  potential  jump  across  the  sheet,  <f>+(t)  +  <j>~{t)  =  2<f>+(t)  — 
A<t>w(r,6,t).  Inserting  this  into  (2.9)  results  in  an  expression  for  the  potential  4>+ : 


4 *<f>+{t)  =  2xA 4>w{ri0,t)+  J  dS  - 

5C,  n*  Uq 

— J  <Iw(t)G(p;q)  dS  +  J  A <j>w(r,6,t) — dS  on  Sew- 

Scw(t)  ScwuSw  q 

(2.10) 


Equations  (2.8)  and  (2.10)  define  the  potential  4>(t)  on  the  blade  surface  beneath 
the  cavity,  Scb ,  and  the  potential  <j>+(t)  on  the  wake  sheet,  5cw(<),  in  terms  of  the 
following  distributions  of  singularities: 

•  continuous  source  and  normal  dipole  distributions  on  SCB  of  strength  g£(t) 
and  </>(t),  respectively 

•  a  source  distribution  on  5cw(t)  of  strength  qw(t) 

•  a  normal  dipole  distribution  on  the  entire  wake  sheet  Sew  U  Sw  of  strength 
A  <j>w{r,9,t). 

On  the  wake  sheet  Sew  U  Sw  the  dipole  strength  A 4>w(r,9,t)  is  converted  along  the 
assumed  wake  surface  with  angular  speed  w,  in  order  to  ensure  that  the  pressure  jump 
in  the  wake  is  equal  to  zero: 


A  <Mr»M) 


A  4>v{r,9,t) 


,  1  »  -  «r(r) 

t<L±£l 


(2.11) 
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where  r,0  are  the  cylindrical  coordinates  of  the  wake  surface  and  0r(r)  is  the  8 
coordinate  of  the  blade  trailing  edge  at  radius  r.  A is  the  unsteady  fully 
wetted  flow  potential  jump  in  the  wake.  A 4>{^(r)  will  be  considered  known;  the 
unsteady  fully  wetted  flow  solution  will  be  found  before  the  unsteady  cavity  solution 
(see  section  4.4). 

Note  that  the  wake  convection  given  by  (2.11)  is  identical  to  the  fully  wetted 
convection  [39].  This  is  allowed  by  the  linearization  of  the  supercavity  in  the  wake, 
where  the  two  force-free  surfaces  are  considered  as  one  (see  equation  (2.5)). 

Everywhere  on  Scb  and  Scw{t )>  either  the  source  distribution  is  known  (a  Neu¬ 
mann  boundary  condition)  or  the  dipole  distribution  is  known  (a  Dirichht  boundary 
condition).  Green’s  formulae  (2.8)  and  (2.10)  may  then  be  discretized  and  rewritten 
as  a  linear  system  of  equations  by  employing  the  boundary  conditions  and  shifting 
all  of  the  known  quantities  to  the  right-hand-side.  This  will  be  described  in  detail  in 
Chapters  3  and  4,  following  the  derivation  of  the  boundary  conditions  in  sections  2.3 
and  2.4. 


2.3  The  Dynamic  Boundary  Condition 

The  dynamic  boundary  condition  (DBC)  requires  that  the  pressure  every  where  inside 
and  on  the  cavity  be  constant  and  equal  to  the  known  cavity  pressure,  pe.  Bernoulli’s 
equation  with  respect  to  the  propeller  fixed  system  is 

7  +  5  l^l*  =  W  +  7  +  51*1*  -  rv  +  9ys-  (2-12) 

where  p  is  the  fluid  density  and  r  is  the  distance  from  the  axis  of  rotation.  Here  qt  is 
the  total  velocity  on  the  cavity  surface.  p0  is  the  pressure  far  upstream  on  the  shaft 
axis;  g  is  the  gravitational  constant  and  ys  is  the  vertical  distance  from  the  horizontal 
plane  through  the  axis  of  rotation,  as  shown  in  Figure  2-1.  ys  is  defined  as  negative 
in  the  direction  of  gravity. 

After  some  manipulation,  and  using  the  definition  of  the  cavitation  number: 


44 


(2.13) 


del-  Po  Pc 

°n  ~  £n2D2  ’ 

where  n  =  ^  and  D  are  the  propeller  revolutions  per  unit  time  and  diameter,  re¬ 
spectively,  we  find  that  the  magnitude  of  the  total  cavity  velocity  satisfies 

|<7t|2  =  n2Z>Vn  [1  -  f(s)}  +  |l/w|2  +  u,2r2  -  2gys  -  2%  (2.14) 

at 

The  function  f(s)  corresponds  to  a  pressure  recovery  law  at  the  trailing  edge  of  the 
cavity  along  the  arc  s  on  the  surface  of  each  spanwise  blade  section.  The  pressure 
recovery  is  given  by  an  algebraic  expression  over  a  specified  portion  of  the  cavity  near 
its  trailing  edge.  This  termination  model  is  described  in  detail  in  Appendix  C. 

Since  the  cavity  boundary  consists  of  two  parts  —  one  coinciding  with  a  portion  of 
the  blade  surface  and  the  other  with  a  portion  of  the  wake  surface  —  the  application 
of  the  dynamic  boundary  condition  on  each  will  be  considered  separately. 

DBC  on  the  Cavitating  Part  of  the  Blade 

In  addition  to  the  expression  (2.14),  the  cavity  velocity  qt  may  also  be  expressed  in 
terms  of  the  directional  derivatives  of  the  perturbation  potential  and  the  components 
of  the  inflow  along  the  same  curvilinear  coordinates  [30].  The  coordinate  system2  on 
the  cavity  surface  consists  of  s  (chordwise)  and  v  (spanwise),  as  shown  in  Figure  2-3: 

„  .  (fc  +  v.)  H  -  d  •  M  +  (g  +  g->  1«  -  <*  •  +  (»  +  U„)  ft  (2,5, 

II*  *  VH2  \9n  / 

with  i  and  v  being  the  unit  vectors  corresponding  to  the  coordinates  s  and  v,  respec¬ 
tively,  and  with  n  being  the  unit  normal  vector  to  the  assumed  cavity.  U„  Uv,  and 
Un  are  the  s,  v  and  n  components  of  the  relative  inflow,  U jn. 

If  s,  v  and  n  were  located  on  the  correct  cavity  surface,  then  the  normal  velocity, 
-f  Un,  would  vanish.  However,  this  is  not  the  case  since  the  cavity  curvilinear 
coordinates  are  approximated  with  those  on  the  propeller  surface.  Nevertheless,  in 

JIn  general  non-orthogonal. 
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applying  the  dynamic  boundary  condition,  the  normal  velocity  is  assumed  to  be 
vanishingly  small.  In  the  fully  nonlinear  scheme,  the  normal  velocity  vanishes  as  the 
solution  converges.  As  shown  in  Chapter  5  and  in  Appendix  B  (for  two  dimensional 
flows),  it  is  found  that  leaving  the  normal  velocity  term  out  of  the  dynamic  boundary 
condition  has  only  a  small  effect  on  the  solution. 

Equations  (2.14)  and  (2.15)  may  then  be  combined  to  form  an  equation  which  is 
quadratic  in  the  unknown  chordwise  perturbation  velocity,  |£.  The  two  solutions  to 
the  quadratic  equation  represent  potential  gradients  which  are  positive  and  negative 
in  the  direction  of  increasing  s.  The  positive  root  is  chosen  to  ensure  that  the  cavity 
velocity  points  in  the  direction  of  increasing  s.  We  can  now  express  in  terms  of 
the  cavitation  number,  the  inflow  velocity,  and  the  unknown  derivatives  |£  and 

~{s,v,t)  =  ~U,  +  cos  9  + 

+  sin0 

(2.16) 

with  0  being  the  angle  between  s  and  v,  as  shown  in  Figure  2-3.  Equation  (2.16)  is 
integrated  once  to  form  a  Dirichlet  boundary  condition  on  <f>: 


The  integral  on  the  right-hand-side  of  (2.17)  is  determined  by  trapezoidal  quadrature, 
as  described  in  section  3.3.  Since  (2.17)  defines  the  strength  of  the  dipole  distribution 
on  the  cavity,  it  may  be  directly  substituted  into  Green’s  formula  2.8. 
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Figure  2-4:  Velocity  diagram  on  the  wake  surface. 


According  to  the  dynamic  boundary  condition  (2.17),  <j>  depends  on  both  its  span- 
wise  and  time  derivatives.  These  terms  will  be  treated  as  knowns  and  will  be  updated 
in  a  time-stepping  scheme  which  will  be  discussed  later.  The  influence  of  the  crossflow 
term  |£  was  studied  first  for  the  case  of  partially  cavitating  3-D  hydrofoils  and  it  was 
found  that  the  global  dependence  of  the  solution  on  the  crossflow  term  was  small. 
This  result  will  be  shown  in  section  4.2.4,  where  it  will  also  be  demonstrated  for  the 
propeller  solution.  The  convergence  of  the  time- derivative  with  iterations,  as  well  as 
its  effect  on  the  solution,  will  be  discussed  in  section  4.4. 

DBC  on  the  Cavitating  Part  of  the  Wake 

The  dynamic  boundary  condition  on  the  cavitating  portion  of  the  wake,  Sew >  may 
also  be  written  as  a  Dirichlet  condition  on  <f>+.  In  this  case,  consider  the  orthogonal 
system  (*,u,n),  shown  in  Figure  2-4.  Assuming  that  a  is  the  direction  of  the  mean 
velocity,  the  total  velocity  on  the  upper  side  of  the  wake  sheet  may  be  written 

V+  =  V+a  +  V+u  +  V+n. 

The  normal  velocity  V+  will  be  omitted  from  the  dynamic  boundary  condition,  as  it 
was  from  (2-4),  with  the  same  justification. 

Applying  Bernoulli’s  equation,  which  is  used  to  define  the  total  velocity  on  the 
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cavity  qt ,  we  have 


(2.18) 


V*  =  Vl««l=  -  (K+)’. 

The  dynamic  boundary  condition  on  Sew  may  thus  be  written 


d*+ 

ds 


-U.  +  ^n*D*<rn  [1  -  f{»)]  4-  \Uw\2  +  -  2 gys  -  2^  -  (V+)».  (2. 


19) 


where  we  have  substituted  the  expression  for  jgt|  given  by  (2.14).  Equation  (2.19) 
may  be  integrated  once  to  form  a  Dirichlet  boundary  condition  on  the  potential  on 
the  upper  wake  surface,  <f>+: 


<f>+(s,u,t)  =  <f>(0,v,t)  +  j 


•TB 


-U.+ 


+^D^<rn  [1  -  f(s)]  +  \UW\2  +  -  2 gys  -  2^  -  (V+)» 


ds.  (2.20) 


The  integral  in  (2.20)  is  also  computed  by  trapezoidal  quadrature.  Equation  (2.20) 
defines  the  potential  <f>+  on  the  upper  side  of  the  wake  and  may  be  directly  substituted 
into  Green’s  formula  (2.10). 

As  with  the  dynamic  boundary  condition  on  the  cavitating  part  of  the  blade,  the 
influence  of  the  crossflow  term,  V^+ ,  was  found  to  be  small.  This  will  be  demonstrated 
in  section  4.2.4. 


2.4  The  Kinematic  Boundary  Condition 

The  kinematic  boundary  condition,  which  is  used  to  determine  the  position  of  the 
actual  cavity  surface  once  the  singularity  strengths  are  known,  is  derived  in  this 
section.  As  in  the  previous  section,  the  cavity  boundary  will  be  divided  into  two 
zones  which  will  be  considered  separately. 

KBC  on  the  Cavitating  Part  of  the  Blade 

The  kinematic  boundary  condition  on  the  cavity  is  the  requirement  that  the  sub- 
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Figure  2-5:  Definition  of  the  cavity  surface  defined  with  respect  to  the  blade  surface, 
stantial  derivative  of  the  cavity  surface  vanishes. 


D 

Dt 


id 

n-h(s,v,t)  =  (—  +  qt  -V)  n-h(s,v,t)  =0  (2.21) 


where  n  is  the  coordinate  normal  to  the  blade  surface  (with  unit  vector  n)  and  h(s,  v,  t) 
is  the  thickness  of  the  cavity  normal  to  the  blade  at  the  point  (s,v)  at  time  t.  Ex¬ 
pressing  the  gradient  in  terms  of  the  local  directional  derivatives 


TT  *)«]  £  +  [*-(*'  *)jl  £  ,  ^  9 

II*  --  +  Udn' 


S  X  V 


(2.22) 


performing  the  dot  product  with  qt  (as  defined  in  (2.15))  and  substituting  the  result 
in  (2.21)  yields  the  following  partial  differential  equation  for  the  cavity  thickness: 


~  [V,  -  cos  9Vv]  +  [Vv  -  cos  9V.)  =  sin3  9(Vn  -  ~) 


(2.23) 


where 


v‘  ~  ds  +  Ul 


v, .  *  d4  +  v, 

ov 


V„ 


«  4 

on 


KBC  on  the  Cavitating  Part  of  the  Wake 

The  kinematic  boundary  condition  on  the  cavity  surface  in  the  wake  may  be 
derived  in  a  similar  fashion,  except  that  now  both  surfaces  of  the  supercavity  must 
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Figure  2-6:  Definition  of  the  cavity  surface  defined  with  respect  to  the  wake  surface. 


be  considered 

^  (n  -9+(s,u,t))  =  J^  +  V+-V  (n-p+(a,u,t))  =0 

£-(n-  g-(3,u,t))  =  -V  (n- <T(s,u,t))  =0  (2.24) 

where  g±(s,u,t)  defines  the  upper  and  lower  cavity  surfaces,  as  shown  in  Figure  2-6. 
Note  that  (s,u,n)  is  an  orthogonal  system.  V+  and  V~  are  the  total  velocities  on 
the  upper  and  lower  sides  of  the  wake  (also  shown  in  Figure  2-6),  respectively,  and 
may  be  written 

V*  =  V*i  +  V±u  +  (2.25) 

The  upper  and  lower  cavity  surfaces,  p(s,  «,<)*,  may  be  written 

g(s,u,t):k  =  C(s,u,t)  ±  ^hw(s,u,t), 

where  C  is  the  cavity  camber  in  the  wake  and  hw  is  the  cavity  thickness.  The  quan¬ 
tities  g ,  C  and  hw  are  all  taken  along  the  normal  to  the  trailing  wake  surface.  These 
quantities  are  also  shown  in  Figure  2-7. 

Expanding  equations  (2.24)  we  find  that 

v+-\™  +  ldM  -  v+  \9c  mJ  +  \dc  i dhj 

n  dt  2  6t  1  da  2  ds  u  2  du 


50 


(2.26) 


dc 

1  dhw 

—  v~ 

ac 

1  dhw 

+  Vu 

dC 

1  dhw 

at 

-  2  ~dt 

—  vt 

ds 

~  2~ds 

du 

2  du 

Taking  the  difference  between  the  two  equations  in  (2.26)  and  assuming  that  s  coin¬ 
cides  with  the  direction  of  the  mean  velocity  so  that 


K+  = 


(2.27) 


yields 


Q  (t)  -  +  2V+  — 

9w{  )  dt  ~  ‘  ds  +  u  du 


(2.28) 


Here  we  have  used  the  definition  of  the  wake  source  strength  (2.7)  and  following 
equalities: 

V*  =  -VJ  and  V*  =  -V'  (2.29) 


which  readily  follow  from  the  assumption  (2.27)  and  the  fact  that  the  free  vorticity 
must  follow  the  mean  velocity  vector. 

To  be  consistent  with  the  dynamic  boundary  condition,  we  assume  that  the  span- 
wise  crossflow  velocity  is  small.  Applying  equation  (2.19),  the  kinematic  boundary 
condition  (2.28)  reduces  to 

[1  -  /(a)]  +  | UW|>  +««r»  -  2gy,  -  2^  =  9.(1)  -  (2.30) 


Note  that  the  cavity  height  on  the  blade  and  in  the  wake,  both  shown  in  Figure 
2-7,  are  defined  differently  and  so  are  given  separate  symbols.  The  computation  of  h 
and  hw  from  (2.23)  and  (2.30)  will  be  described  in  section  4.1.3. 

The  position  of  the  cavity  surface  over  the  blade  surface  is  determined  by  adding 
the  cavity  thickness  h  normal  to  the  blade  surface  at  the  midspan  of  the  panel  bound¬ 
aries.  In  the  wake,  the  cavity  camber  C(s,u,<)  must  first  be  determined.  An  expres¬ 
sion  for  C(s,u,t)  may  be  found  by  adding  the  two  equations  (2.26)  and  dividing 
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Figure  2-7:  Definition  of  the  cavity  camber  and  height  for  a  supercavitating  section 
of  the  propeller  blade. 

through  by  two: 


(2.31) 


where  Un  is  the  inflow  velocity  normal  to  the  wake  sheet.  Equation  (2.31)  is  numeri¬ 
cally  integrated  to  determine  the  camber  surface  in  the  wake.  At  the  trailing  edge  of 
the  blade,  the  continuity  of  camber  and  thickness  is  imposed: 


<.  = 


(2.32) 


Here,  the  superscripts  +  and  —  denote  just  upstream  and  just  downstream  of  the 
trailing  edge,  respectively.  C£  is  the  value  of  the  camber  just  upstream  of  the  trailing 
edge.  It  is  determined  by  adding  |  to  the  trailing  edge  along  the  blade  normal  (see 
Figure  2-7).  The  quantity  h~u  is  determined  by  interpolating  the  upper  cavity  surface 
over  the  blade  at  the  trailing  edge  and  computing  its  normal  offset  from  the  wake 
sheet.  The  upper  and  lower  surfaces  of  the  cavity  in  the  wake  are  then  determined  by 
adding  and  subtracting  half  of  the  cavity  thickness  hw  from  the  camber  surface.  This 
defines  the  cavity  surface  at  the  midspans  of  all  the  spanwise  strips.  The  surface  of 
the  cavity  at  the  strip  boundaries  are  determined  by  interpolation  and  extrapolation. 


KBC  on  the  Wetted  Part  of  the  Blade 


The  kinematic  boundary  condition  on  the  wetted  portion  of  the  blade,  Sws ,  defines 
the  source  strength  there  in  terms  of  the  known  inflow  velocity: 

Q  J 

(0  =  —  ^in(a'9’  y**  9  €  5ws(0  (2.33) 

where  xq,yq,zq  are  the  coordinates  of  the  point  with  respect  to  the  propeller  fixed 
system.  As  in  the  case  of  fully-wetted  flow,  this  boundary  condition  may  be  directly 
substituted  in  Green’s  formula. 


2.5  The  Kutta  Condition 

The  Kutta  condition  requires  that  the  fluid  velocity  be  finite  at  the  blade  trailing  edge. 
It  was  found  necessary,  in  the  fully  wetted  steady  flow  propeller  solution,  to  satisfy  a 
nonlinear  Kutta  condition  which  ensures  pressure  equality  between  the  suction  and 
pressure  sides  at  the  trailing  edge  [47].  The  so-called  pressure  Kutta  condition  was 
later  refined  by  Kinnas  and  Hsin  [38],  who  developed  an  efficient  iterative  solution. 
In  the  present  work,  an  extension  of  Morino’s  steady  Kutta  condition  is  applied. 
Development  of  a  pressure  Kutta  condition  for  the  cavitating  propeller  is  left  for  the 
future. 

The  value  of  the  dipole  strength,  A <£r(r,f),  at  the  blade  trailing  edge  at  time  t  is 

A <h{r,t)  =  #r(r,t)  -  ^(r.t)  =  T(r,t)  (2.34) 

where  4>t  an<^  4>t  are  the  values  of  the  potential  at  the  blade  trailing  edge  at  radius 
r  on  the  suction  side  and  the  pressure  side,  respectively.  The  potential  jump  there  is 
also  equal  to  the  circulation  T  at  time  t  around  the  blade  section  at  radius  r.  This 
condition  is  equivalent  to  requiring  the  shed  vorticity  from  the  blade  trailing  edge 
to  be  proportional  to  the  time  rate  of  change  of  the  circulation  around  the  blade 
(Kelvin’s  law).  An  extension  of  Morino’s  Kutta  condition  for  steady  flow  [52],  in 
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which  the  potential  jump  at  the  trailing  edge  of  the  blade  is  simply  replaced  by  the 
potential  jump  at  the  nearest  control  points,  is  applied.  This  is  addressed  in  detail 
in  section  4.4. 
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Chapter  3 


Discrete  Formulation 


The  greatest  complexity  of  the  cavity  flow  problem  is  the  fact  that  the  cavity  extent 
(planform)  is  not  known  a  priori ;  its  computation,  in  fact,  is  a  primary  goal  in  the 
solution.  Clearly,  since  the  application  of  the  boundary  conditions  relies  on  knowing 
the  extent  of  the  cavity,  computing  it  will  require  some  sort  of  iterative  process  of 
successive  guesses.  In  this  chapter,  a  brief  introduction  to  the  iterative  procedure  will 
be  given,  followed  by  a  detailed  description  of  the  discretized  integral  equations  and 
boundary  conditions. 


3.1  Solution  Methodology 

The  local  cavity  length  (defined  as  the  arclength  of  the  projection  of  the  cavity  on  the 
nose-tail  helix)  is  given  at  each  radius  r  by  the  function  l{r,t).  For  a  given  cavitation 
number,  <rn,  the  cavity  planform  /(r,  t)  will  be  determined  from  the  requirement: 

$( J(r,  t),  r;  <rn)  h(l(r,  t),  r,  <)  [or  h»(/(r,  t ),  r,  <)]  =  0.  (3.1) 

Equation  (3.1)  requires  that  the  cavity  close  at  its  trailing  edge;  this  requirement  will 
be  used  as  the  basis  of  an  iterative  solution  to  find  the  cavity  planform. 

The  iterative  method  for  finding  the  cavity  length  for  given  cavitation  number 
evolved  from  the  numerical  solution  of  the  fixed-length  twn.dimen«on»l  flow  problem, 
in  which  the  cavity  length  on  a  two-dimensional  hydrofoil  is  assumed  to  be  known 
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while  the  cavitation  number  is  (an  additional)  unknown  [35].  In  that  case,  a  cavity 
closure  condition  —  requiring  that  the  cavity  close  at  its  trailing  edge  — is  imposed 
to  provide  the  necessary  additional  equation.  However,  when  solving  the  so-called 
“direct”  problem,  where  the  cavitation  number  is  known  and  the  cavity  length  is 
to  be  determined,  there  is  one  fewer  unknown  and  thus  the  closure  condition  may 
be  discarded.  The  resulting  indeterminacy  is  then  experienced  in  the  following  way: 
solving  the  boundary  value  problem  for  given  cavitation  number  by  guessing  the 
extent  of  the  cavity  results  in  a  cavity  which  does  not  close  unless  the  cavity  length  is 
the  one  which  corresponds  to  the  given  cavitation  number.  Put  another  way,  for  given 
cavitation  number  the  correct  cavity  length  is  the  one  for  which  the  cavity  closes. 
This  reveals  the  foundation  of  the  algorithm  for  finding  the  correct  cavity  length:  the 
solution  is  found  for  successive  guesses  of  the  cavity  length  until  the  solution  finds  a 
closed  cavity1.  Note  that  the  Newton-Raphson  iterative  method  for  root-finding  is 
ideal  for  this  application. 

Figure  3-1  depicts  the  process  of  finding  the  correct  cavity  length  for  a  partially 
cavitating  2-D  hydrofoil.  Shown  are  a  succession  of  cavity  solutions  for  a  NACA16006 
hydrofoil  at  an  angle  of  attack  a  =  4°  corresponding  to  a  succession  of  guesses  of  the 
cavity  length.  The  cavitation  number  for  all  solutions  is  <r  =  1.097*.  Notice  that 
the  cavity  closes  for  the  guess  l  =  0.4,  as  expected.  Also  shown  in  Figure  3-1  is  the 
openness  of  the  cavity  at  its  trailing  edge,  6,  plotted  against  the  cavity  length,  l.  The 
solution  of  the  direct  problem  amounts  to  finding  the  zero-crossing  of  this  curve. 

The  iterative  solution  for  the  cavity  extent  may  also  be  applied  in  three  dimen¬ 
sions,  where  several  difficulties  are  encountered.  The  most  obvious  difficulty  is  the 
need  to  determine  the  cavity  length  in  both  chordwise  and  spanwise  directions.  An 
additional  challenge  is  the  accurate  prediction  of  the  cavity  planform  without  regrid- 
ding  the  blade  surface  for  each  iteration (*.  Perhaps  the  most  severe  difficulty  arises 

1  An  open  cavity  model  is  also  very  easily  implemented  with  this  method.  Instead  of  requiring 
the  cavity  to  dose,  the  desired  thickness  of  the  cavity  at  its  trailing  edge  may  be  specified. 

•This  cavitation  number  corresponds  to  the  fixed-length  solution  for  l  =  0.4. 

•The  solution  for  each  cavity  length  shown  in  Figure  3-1  was  found  by  regridding  the  hydrofoil 
for  each  guess  of  the  cavity  length  so  that  the  cavity  always  ended  at  a  panel  boundary. 
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Figure  3-1:  NACA16006  hydrofoil  at  a  =  4°  and  <r  =  1.097  (corresponding  to  the 
correct  length  l  =  0.4).  Left  plot  shows  the  predicted  cavity  shapes  for  different 
guesses  of  l.  Right  plot  shows  the  resulting  behavior  of  6  vs  l. 

from  the  need  to  simultaneously  treat  the  occurence  of  both  partial  and  supercavita¬ 
tion.  Before  providing  answers  to  these  difficulties,  however,  it  is  useful  to  describe 
the  solution  for  a  given  guess  of  the  cavity  planform.  In  the  following  sections,  I  will 
describe  in  detail  the  numerical  discretization  of  the  boundary  value  problem  which 
has  been  detailed  in  Chapter  2. 

3.2  Discretized  Green’s  Formula 

The  objective  of  the  numerical  analysis  is  to  invert  equations  (2.8)  and  (2.10)  subject 
to  the  kinematic  boundary  condition  (2.33),  the  dynamic  boundary  conditions  (2.17) 
and  (2.19),  and  the  Kutta  condition.  To  accomplish  this,  we  first  tag  one  blade 
with  the  label  key  blade.  The  solution  at  a  given  time  step  will  be  obtained  only 
for  the  key  blade,  with  the  influence  of  the  other  blades  corresponding  to  an  earlier 
solution  of  the  key  blade.  This  time  stepping  algorithm  will  be  detailed  in  section  4.4. 
The  key  blade  surface  is  discretized  into  N  chordwise  and  M  spanwise  quadrilateral 
panels  with  the  corners  lying  on  the  blade  surface  Scb  and  with  the  control  points 
located  at  the  panel  centroids.  An  example  of  a  discretized  blade  is  shown  in  Figure 
3-2.  The  source  and  normal  dipole  distributions  on  each  panel  are  approximated 
with  constant  strength  distributions.  The  trailing  wake  is  discretized  into  panels  at 
constant  angular  intervals  =  u>A t  with  At  being  the  time  step.  The  blade  and 
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Figure  3-2:  Discretization  of  the  propeller  blade,  the  cavity,  and  their  trailing  wakes; 
N=16,  M=9. 

trailing  wake  discretization  is  identical  to  that  in  the  case  of  folly  wetted  unsteady 
flow  [22,  39,  38].  If  we  call 

Nws  =  Number  of  wetted  panels  on  the  blade 
Ncb  =  Number  of  cavitating  panels  on  the  blade 
New  =  Number  of  cavitating  panels  in  the  wake, 

then,  among  the  discrete  sources  and  dipoles,  we  have 

•  Nws  known  source  strengths,  via  (2.33) 

•  Ncb  known  dipole  strengths,  via  (2.17) 

•  New  known  dipole  strengths,  via  (2.20). 

The  following  are  then  unknown  and  must  be  solved  for: 

•  Nws  dipole  strengths  on  the  wetted  blade  surface. 

•  Ncb  source  strengths  on  the  cavitating  blade  surface,  and 

•  New  cavity  source  strengths  in  the  super  cavitating  wake. 
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Prior  to  substituting  the  expressions  for  the  known  singularity  strengths,  the  dis¬ 
cretized  Green’s  formulae  (2.8)  and  (2.10)  appear  as  follows: 


r*M  M  (  N  r  suh 

IE  E  U»A.(«>  -  BL.%  (ft) 

fcsl  mil  ln=l  L  OTlnm 

Nw 
Iml 


cJm«L(ft)+ 


i  =  1,  ...,1V  x  M  (3.2) 


and 


4x0*(n)  = 


E 


l=i 


fLqU*)+ 


Nw  h  ) 

+  EW3mA<Jft)|  *  =  1 . New 


(3.3) 


where  n  is  the  discrete  time  step,  Nw  is  the  number  of  panels  on  each  strip  of  the 
wake,  and  NCm  is  the  number  of  cavitating  panels  on  each  strip  of  the  wake. 

In  equations  (3.2)  and  (3.3),  the  indices  t  and  j  map  quantities  to  panels.  For 
example,  &(n)  is  the  potential  at  the  control  point  on  the  ith  panel  at  the  nth  time 
step.  However,  each  panel  may  also  be  identified  as  the  nth  panel  on  the  strip, 
so  <t>i  may  also  be  written  In  what  follows,  these  indexing  alternatives  will  be 
used  interchangeably  to  maximize  compactness. 

Note  that  the  subscript  q  has  been  dropped  from  references  to  the  discretized 
potential,  and  the  source  strength,  since  these  quantities  are  defined  at 

the  control  points  which  are  by  definition  variable  points  in  the  integration.  For  the 
remainder  of  this  section,  all  of  the  flow  variables  will  be  assumed  to  be  defined  at 
the  current  time  step  and  n  will  be  dropped  from  the  equations. 

and  Bj^  are  the  potentials  induced  at  the  itH  control  point  on  the  key  blade 
by  a  unit  strength  dipole  and  a  unit  strength  source  at  the  nth  panel  on  the  m*  strip 
of  the  k**  blade  (k  =  1  refers  to  the  key  blade).  When  the  Ith  control  point  lies  on 
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the  ntk  panel  on  the  mth  strip  of  the  key  blade,  then 

=  2ir  fori. — >(n,m). 

is  the  potential  induced  at  the  itk  panel  on  the  key  blade  due  to  a  unit  strength 
source  at  the  nth  panel  on  the  mth  strip  of  the  wake  of  the  ktk  blade.  is  the 

potential  induced  at  the  itk  control  point  on  the  key  blade  by  a  unit  strength  dipole 
at  the  Ith  panel  of  the  mth  strip  of  the  wake  of  the  kth  blade.  Similarly,  D*nm  and  Eknm 
are  the  potentials  induced  at  the  itk  wake  control  point  by  a  unit  strength  dipole  and 
a  unit  strength  source  at  the  nik  panel  of  the  mtk  strip  of  the  kih  blade.  Fknm  is  the 
potential  induced  at  the  iik  wake  control  point  due  to  a  unit  strength  source  at  the 
ntk  wake  panel  on  the  mtk  wake  strip  of  the  kth  blade.  is  the  potential  induced 
at  the  ith  wake  panel  by  a  unit  strength  dipole  lying  on  the  mth  strip  of  the  wake  of 
the  kth  blade,  where 

W?nm  =  2*  for  *  < — ♦  (n,m). 

The  shape  of  the  surface  bounded  by  the  edges  of  each  quadrilateral  panel  is  ap¬ 
proximated  by  a  hyperboloidal  surface,  and  the  corresponding  influence  coefficients 
are  determined  analytically.  The  need  for  hyperboloidal  panels  was  found  to  be  nec¬ 
essary  for  the  convergence  and  consistency  of  the  steady  flow  propeller  solution,  espe¬ 
cially  when  applied  to  extreme  geometries  [38,  53,  21].  Discussion  of  the  computation 
of  these  influence  coefficients  may  be  found  in  [55]  and  [22]. 

Equations  (3.2)  and  (3.3)  may  be  rewritten  to  reflect  the  fact  that  at  any  given 
time  step  only  the  potentials  on  the  key  blade  are  unknown,  while  the  rest  are  assumed 
to  be  known.  Equation  (3.2)  becomes 

E  l  E  knm* nm  -  Binm^  ]  -  £  CinmQnm  +  1  =  {rHs} 

mxl  [rml  ^  *nn”*J  n=l  J  1  >  SOBi 

i  =  1,...,  JV  x  M  (3.4) 
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Figure  3-3:  Indexing  system  for  the  discrete  equations. 


18 


d<p 

dnjm 


n  =  1, Nle ;  Nle  +  NcAVm, ...» N 
m  —  1, M 


(3.8) 


where  U-n^  is  the  inflow  velocity  at  the  current  time  step,  defined  at  the  n**  control 
point  on  the  mth  spanwise  strip  of  the  blade.  The  system  of  indexing  is  shown  in 
Figure  3-3. 

A  single  discrete  equation  may  be  written  to  replace  the  dynamic  boundary  con- 


and  (3.3)  becomes 


=  4  *4>i  + 


i  —  1*  •••)  New  (3.5) 


where 


__  f>k 
“inm 


d£h 

dtlnm 


**em 


_  V'  c*  o* 

/  ,  ^inmx  r 


n=l 


Nw 

+  E«4i^ 

1—2 


»lm 


! 


(3.6) 


and 


[rhs} 

■  '  Scwi 


ry h  jjk  __  pk  ^ 

^inmTnm  ciiim  Qrinr 


*©, 

n=l 


mm^nm 


Nw 

+ T. "'*># 

1-2 


h 

toim 


)• 


(3.7) 


In  (3.4)  and  (3.5)  the  superscript  k  =  1  is  implied.  In  (3.6)  and  (3.7)  the  potential 
^nm  i*  «qual  to  <f>hm  taken  from  the  key  blade  solution  at  a  previous  time  step,  cor¬ 
responding  to  the  current  location  of  blade  k.  The  same  equivalence  is  true  for  the 
source  strengths  and  Q*m  and  the  wake  dipole  strengths 


3.3  Discretized  Boundary  Conditions 

The  boundary  conditions  may  now  be  discretized  and  incorporated  in  Green’s  formu¬ 
lae  (3.4)  and  (3.5). 

The  discretized  version  of  the  kinematic  boundary  condition  on  the  blade  (2.33) 
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ditions  on  the  blade  (2.17)  and  in  the  wake  (2.20) 


4>(NLM+n)m  =  <t>0m  +  ^nm 


n  =  l,...,i\TCilVm  -f  NCm 
m  =  1, ...,  M 


(3.9) 


where 


$  4*  / 

« nm  —  i 


T  [If  fc°m  (2*16)]  da  U~  1*  •••»  NcAVm 

0  L  J 

/  [§£  from  (2.16)]  da  +  /  [§£  from  (2.19)]  da 


o 

•T  9m 


<nm 

/ 

•TBm 

n  =  Wcmv*  +  1,  — ,  NCAVm  +  ^c„ 


(3.10) 


m  =  1, M 


and  ste„.  is  the  value  of  a  at  the  trailing  edge  of  the  blade  on  the  mth  span  wise  strip, 
as  shown  in  Figure  3-4.  The  integrals  in  equation  (3.10)  are  computed  by  trapezoidal 
quadrature.  The  values  of  the  integrands  are  computed  at  the  control  points  and 
at  the  leading  edge  of  cavity,  where  a  =  0.  The  arclength  between  two  consecutive 
control  points  is  approximated  by  the  sum  of  the  linear  distances  between  the  control 
points  and  the  panel  edges  (see  Figure  3-4) 

1  n_1  i 

anm  —  o(ASjm  +  Aj(j+i)m) 

1  ;=i 1 


The  quantity  4>q  in  (3.9)  is  the  perturbation  potential  at  the  detachment  point  of 
the  cavity.  It  is  also  an  unknown.  This  term  is  expressed  as  a  cubic  extrapolation  in 
terms  of  the  unknown  potentials  on  the  wetted  panels  on  the  same  strip  adjacent  to 
the  cavity  detachment.  This  is  shown  schematically  in  Figure  3-5. 

Locally,  the  potential  is  assumed  to  fit  a  3r<*-order  polynomial  in  arc-length  of  the 
form 

<f>{a)  —  Aa 8  +  Ba 2  +  Ca  +  D  (3.11) 

where  the  coefficients  A,  £,  C,  D  may  be  written  in  terms  of  the  local  unknown  po- 
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Figure  3-4:  Definition  of  the  arclength  on  a  sp&nwise  strip,  used  in  the  trapezoidal 
integration  of 

tentials  and  the  cavity  velocity  by  applying  the  following  boundary  conditions  (see 
Figure  3-5): 

.  4(1  =  0)  ^  4>z„ 

•  <f>(3  =  S2)  d=  <t>2m 

•  4(3  =  Si)  =f  <t>in 

•  f»(^  =  to)  =  qom  =f  §£  at  the  cavity  leading  edge  (from  2.16). 

Solving  for  the  coefficients  in  (3.11)  in  terms  of  these  quantities  yields  the  following 
definition  of  <f>om- 

<t>Om  =  +  Rimfam  +  ^3m<^3m  +  R*m  qO*.  (3.12) 

where  Rim, ...,  are  rational  functions  of  So,  3 1,  and  s2: 

Rlm  =  —  [1  —  25o*i)] 

jl 

R}m  =  — ”  [2iori$r  —  23otiT2  +  r2] 

7*1 

3 a 

Rzm  =  1  +  —  [2ritjio  4-  2<i5o(l  —  r*)  +  t*j  —  1] 

ri 

R*m  =  [— (<i  —  ^2rz)  +  Sjtjl - ”(^i  —  Wi)  + 

iri  J  ri 
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Figure  3-5:  Schematic  representation  of  the  treatment  of  No¬ 


where 


rj  =  £2^(12  -  2So)  —  +  2J0si 

-  £l  **  ~~ 

Tj  S2  H  -  35o 

,  _  ~  2J° 

1  S§  - 

t,  =  * 

5|  —  3jJs2 


4>o  is  thus  defined  in  terms  of  local  potentials,  which  are  unknown,  and  the  cavity 
velocity,  which  is  known.  As  a  result,  <f>o  is  split  so  that  the  known  part  is  on  the 
right-hand-side  and  the  unknown  part  is  on  the  left-hand-side. 

For  a  discussion  of  the  importance  of  to  the  solution,  see  in  section  4.3. 

Incorporating  the  discrete  boundary  conditions,  the  discretized  Green’s  formulae 
take  the  following  form: 


On  Scb 


M  (  Nlb  ^ 

£  |  £  Ainmtnm  +  £  A^m4>nm  +  +  -Rjm&n.  + 

m=l  \n=l  nsWtB+WcjiVm+l 


I*LM+NcAVm  QQ 


£  -  £  CiimC?/m  +  Wlm  [(1  -  SOPm)^m 

_^1_  _  c/Wtifn  * _ . 

WLJS4-1  ‘=1 


n=NX,£+l 


-  — 

M  (  NCAVm  1*LB  V'  D  rr  A 

=  V  |  -  £  ^(AriJI+n)m$nm  ~£  BinmUmnm  -nnm-  £  «inm ^  •  ««•» 
^i(  n=l  n=l  tisNj,g  +1 

-  SOPmWiim«Nm  -  «4.qom(im|  +  {RHS}ScJ|. 


*  =  l,...,JVxM  (3.13) 


where 


&m  — 


JVjta+Novm 

j4inm  +  SOPmWiim. 

n=Nj,a+l 


On  Sew 


M  (  Nlb  ^  i  u 

£  {  £  Dinm<t>nr*  +  £  Dinm^nm  + 

m=l  (  n=l  ns=i^i,*+JVcitVm+1 

T **»»  -  E  »„<?.-  +  »s-  [d  -  sop.»*.  -*.]}- 
„Jr.+.  9n“"-  >-■  J 


M  f  *CAVm  *£* 

-4*¥i  +  £  {  -  £  Di(SLB^)m^nm  “£  *  *W 

m=l  V  n=1  n=1 


iV 


£  EinwVin*m  ' 
n*Wt,a+tfc4Vm+1 


-S0Pmw?am«wm  -  «4«qo«.u|  +  {RHS}Scwr4 


*  =  1, ...,  New  (3.14) 


where 
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NLM+NcAVn 

53  Anm  +  S0PmW7lm  4-  4tt  if  t  is  on  the  mth  strip 

nsiST^j+l 
Nlb+Ncav,* 

53  A«m  +  SOPmW?lm  otherwise 

n=Nia+l 

In  order  to  write  (3.13)  and  (3.14),  which  are  true  for  arbitrary  cavity  planforms 
including  partial  and/or  supercavitation,  it  was  necessary  to  introduce  the  switch 
SOPm  to  indicate  whether  a  given  strip  is  supercavitating  or  partially  cavitating. 
This  switch  is  defined  as  follows: 


SOP 


1.0  if  the  mth  strip  is  supercavitating 
0.0  otherwise. 


Chapter  4 


Solution 


The  aim  of  this  chapter  is  to  complete  the  description  of  the  solution  methodol¬ 
ogy.  The  system  of  equations  which  must  be  solved  for  each  planform  iteration  at 
each  time  step  is  written  down  and  the  efficient  and  robust  algorithm  for  arriving  at 
the  correct  planform  is  described.  Issues  attendant  to  the  algorithm  are  discussed, 
including  uniqueness  and  sensitivity  to  tolerances.  The  method  of  time-marching 
is  outlined,  followed  by  the  details  of  obtaining  post-processed  results,  such  as  the 
pressure  distribution  and  the  cavity  shape. 


4.1  Solution  for  Given  Cavity  Planform 

As  mentioned,  the  planform  of  the  cavity,  /(r,£),  is  unknown  and  has  to  be  deter¬ 
mined  as  a  part  of  the  solution.  However,  since  l(r,  t)  will  be  determined  iteratively 
by  successive  guesses,  the  problem  of  determining  the  cavity  shape  at  a  prescribed 
(“guessed”)  planform  will  be  addressed  first. 

For  a  given  guess  of  the  cavity  planform,  (3.13)  and  (3.14)  are  satisfied  at  the 
control  points  on  the  blade  and  in  the  wake  and  the  resulting  linear  system  is  solved 
with  respect  to  the  unknown  source  and  dipole  strengths.  However,  before  solving  the 
system  of  equations,  I  have  to  address  the  following  problem:  for  a  given  discretiza¬ 
tion  and  an  arbitrary  guess  of  I(r,  t),  the  trailing  edge  of  the  cavity  will  not  necessarily 
coincide  with  a  panel  boundary  in  the  chordwise  direction.  This  is  a  problem  because 
the  singularity  distributions  which  represent  the  cavity  span  an  integral  number  of 
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panels.  There  are  several  options  for  circumventing  this  problem,  one  being  to  simply 
re-panel  the  blade  surface  so  that  the  cavity  trailing  edge  does  coincide  with  a  panel 
boundary.  This  idea  is  rejected  immediately  due  to  the  inordinate  expense  of  recom¬ 
puting  influence  coefficients1.  Another  option  is  to  approximate  the  cavity  planform 
by  ending  the  cavity  at  the  nearest  panel  boundary  on  each  spanwise  strip.  Although 
this  method  has  been  applied  with  some  success,  it  was  found  that  the  convergence  to 
the  correct  cavity  planform  suffered  due  to  discontinuous  intermediate  discrete  cav¬ 
ity  planforms  [36].  To  address  this  problem  we  invented  the  split-panel  technique,  in 
which  the  panel  at  the  trailing  edge  of  the  cavity  is  split  into  a  cavitating  and  a  non- 
cavitating  part.  Figure  4-1  describes  this  schematically.  In  the  figure,  the  discretized 
blade  is  represented  by  a  rectangular  grid  and  the  smooth  cavity  planform  overlays 
the  grid  and  is  seen  to  cut  one  chordwise  panel  on  each  spanwise  strip.  The  lengths 
of  the  split  panels  at  their  midspans  are  lim  and  llU,  respectively.  The  split-panel 
method  allows  the  discrete  cavity  planforms  to  be  smooth  without  the  added  burden 
of  re-computing  influence  coefficients. 


4.1.1  Details  of  the  Split-Panel  Method 

The  source  and  dipole  strengths  on  the  split  panels  ( <f> ^  and  §t  )  on  the  blade 
are  defined  as  averages 

=  thshsLttSslSs.  on  the  split  nancl 


llm  -I-  Itu 


J  '  _  tn\L,n  onnm.  _ it. _ l.*i _ 1 

= - - — ; -  Q-  = - : — —7 -  on  the  split  panel 

*£«+*«*  iLm+tRm 

(4.1) 

with  the  weights  lin  and  Ijtm  taken  as  the  lengths  at  the  midspan  of  the  split  pan¬ 
els.  For  compactness,  the  split-panel  dipole  and  source  strengths  will  henceforth  be 


written 


= + *L = aLFi~ + tLF*-  (4-2) 


1  Remember  the  theme  discussed  in  section  2.2:  avoid  re-computing  influence  coefficients! 
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Figu«  44:  Schematic  of  a  cavity  planform,  ehowmg  the  .pfi,  pueU. 


mini 

Figure  4-2:  Sketch  defining  split  panel  extrapolations, 
where  the  factors  Fim  and  Fr^  are  defined  as  follows: 


FLm 


iRm 
+  *  JU 


The  quantities  <pim  and  are  known  from  the  boundary  conditions,  (2.17) 
and  (2.33).  The  quantities  (pR^  and  are  unknown,  but  may  be  defined  as 
extrapolations  of  the  unknowns  on  adjacent  panels,  similar  to  the  definition  of  <£om  in 
section  3.3.  Using  the  notation  of  Figure  4-2,  <Pr^  is  defined  by  a  cubic  extrapolation 
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of  the  unknown  potentials  at  the  four  adjacent  control  points: 


4>R~  =  P\m4>\m  +  Pimfam  +  +  P (4.3) 


where 


Pu 

P*~ 

P^ 


S|  ~  Sg(g3  +  5a)  +  3oSjS8 

h(h  -si)(*3  ~  *i) 

-Sg  -  3g(5i  -f  s3)  -  spii S3 
Sj(S$  —  ij)(Sj  —  Si) 

Sp  ~  H(h  +  Si)  -f  SoSii2 

S3(S3  —  S2)(S3  —  Si) 

j  _l_  ~Sq  ~f~  Sq(Si  +  S3  +  h)  —  So(3iia  +  jii,  4-  s3s3) 

SiSjS) 


(4.4) 


|^L  is  also  defined  as  a  polynomial  extrapolation,  but  I  found  it  necessary  to  include 
a  term  which  is  square-root  singular  at  the  cavity  trailing  edge.  In  this  way,  the  correct 
behavior  of  |£(s)  near  the  end  of  the  cavity  is  captured: 


3j~(S)  2;  a(i  -  3L)2  +  b(S  -  Si,)  +  +  d. 

On  y/a  —  si 

This  allows  me  to  write  ^Lm  as  (again  using  the  notation  of  Figure  4-2) 

d<j>  „  d<f>  „  d<f>  _  d<f>  _  d<ft 
=  A,m8n.m  +  +  + 


(4.5) 


where 


S  ji  f»iTi  +  SjTi  Qi  »,1  ffcT1+T, /(5.T,  s 

D  t7  tTJ  ^  VtT  Q,l 

■ 

SiTi  -f  S$Ta 

Dy/Tl 

o  -»  /  Sj  SiTi-f^Tj^  _  (  1 

5l-  “  (tT  + - D - )  ~  n  (~h 


,  S»T3  S3T3QiT3 
S,Ti  isTiD 


S»TaQ3\ 
ssD  ) 
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is  —  h  .  NQi 


1  ,  (h-h)Ts  .  NQ1T3  NQs 


and  where 


<3  —  <1 

y/it  —  *3  V*l  —  is 
O  —  _  3a  3}  —  <i 

y/*L  -  5a  Vii  -  *i  Vft 

Q*  =  IP  ~  ~jL 

y/SL  -  *3 

Ti  =  3 j3j  -  3j3S 
Tj  =  i\i\  —  SjSj 

T,  =  3^-2313, 

N  =  (ii  -  ia)Ta  +  (h  -  ia)T, 

D  =  QjTj-QiT,  (4.7) 

In  a  similar  manner,  the  source  strengths  at  the  wake  panels  which  are  split  by 
the  trailing  edge  of  the  supercavity  are  defined  as 


=  &%■  <“> 

Qim  is  defined  as  a  polynomial  extrapolation  in  terms  of  the  adjacent  unknown  cavity 
source  strengths,  similar  to  the  definition  of 

Qi*m  —  +  Timqlm  +  Tsmtom  +  Tlbnft*.*  (4.9) 


The  expressions  which  define  7\m, ...,  r4m  are  identical  to  the  expressions  which  define 
5im, ...,  S4m  in  (4.6).  However,  since  3  for  the  Qim  extrapolation  is  different  from  the 
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5  used  in  the  t  extrapolation,  separate  symbols  are  necessary. 

Error  Due  to  Split-Panel  Method 

There  are  two  sources  of  error  in  the  split-panel  method.  They  are: 

•  due  to  the  fact  that  4>Rm  and  are  defined  by  extrapolations  which  are  not 
exact 

•  due  to  the  averaging  of  the  influence  of  two  panels  into  one. 

The  first  source  of  error  is  bounded  by  that  of  the  extrapolation,  since  the  behav¬ 
ior  of  the  source  strength  is  square-root  singular  near  the  cavity  trailing  edge  (while 
the  dipole  strength  is  smooth).  The  correctness  of  the  extrapolation  may  be  gauged  in 
the  following  way:  first  solve  for  the  singularity  strengths  for  a  two-dimensional  flow 
for  given  <r,  a  given  guess  of  the  cavity  length  (without  searching  for  the  closed  cav¬ 
ity)  and  a  given  discretization  (using  the  split-panel  method).  Next,  solve  the  same 
problem  where  the  split  panel  is  treated  as  two  distinct  panels.  The  source  strength 
on  the  left  split  panel,  |££,  should  be  close  to  the  computed  source  strength.  In  ad¬ 
dition,  the  source  strengths  on  neighboring  panels  from  the  two  computations  should 
be  equivalent.  Figure  4-3  shows  the  two  computations,  where  the  is  defined  by 
a  plain  cubic  extrapolation  of  neighboring  unknowns  (left)  and  by  the  extrapolation 
just  described  (right).  Comparing  the  two  plots,  it  is  clear  that  the  extrapolation 
with  the  square-root  singular  term  is  necessary. 

The  second  source  of  error,  due  to  the  averaging  of  the  panel  influence,  is  inves¬ 
tigated  by  the  following  test:  consider  the  influence  of  a  unit  source  panel  in  two 
dimensions  which  has  been  split  into  two  panels  at  some  arbitrary  fraction  of  its 
length.  As  depicted  in  Figure  4-4,  the  total  influence  at  some  field  point  is  the  sum 
of  the  influence  of  the  left  and  right  interior  panels.  Denoting  Bw  as  the  influence  of 
the  original  unsplit  panel,  the  induced  potential  B  may  be  written 

*-*,+*-*[£+ f£] 
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Figure  4-3:  Comparison  of  the  source  strength  computed  with  a  cavity-fit  grid  to 
that  computed  by  the  split  panel  method.  In  the  plot  on  the  left,  |££  is  defined  as  a 
cubic  extrapolation  of  neighboring  unknowns.  On  the  right,  |££  is  defined  as  a  cubic 
extrapolation  with  singular  term,  as  in  the  text.  Computations  are  for  a  NACA16006 
foil  with  <r  =  1.0,  a  =  3°,  and  a  guess  of  the  cavity  length  f  =  0.4. 

where  BL  and  BR  are  the  induced  potentials  due  to  the  left  and  right  interior  panels, 
respectively.  Referring  to  equation  (4.2),  it  is  evident  that  the  split-panel  method 
amounts  to  approximating  the  quantities  and  -jfe  with  the  constants  F&  and  Fr. 
In  other  words,  for  all  field  points  the  split-panel  method  makes  the  approximations 
-J^  fa  Fl  and  «  Fr.  This  approximation  may  be  expected  to  do  well  if  the  field 
point  is  far  from  the  singularities,  but  break  down  for  field  points  nearby.  In  Figure 
4-4,  the  ratio  ■§£  is  shown  for  a  series  of  field  points  arrayed  along  a  line  adjacent 
and  parallel  to  a  split  source  panel.  To  be  consistent  with  the  hydrofoil  solution,  the 
closest  field  point  is  one  half  a  panel  length  away.  From  this  figure,  it  is  clear  that 
the  approximation  breaks  down  within  several  panel  lengths. 

The  effect  of  the  split  panel  error  on  the  global  solution  may  be  seen  by  comparing 
the  cavity  length  for  fixed  geometry  and  a  given  <r  obtained  from  various  initial 
guesses.  If  the  split  panel  error  has  a  large  effect  on  the  solution,  then  the  final  cavity 
length  (corresponding  to  a  closed  cavity  at  the  given  cavitation  number)  would  depend 
strongly  on  the  initial  guess.  It  has  been  found  that  the  cavity  length  depends  very 
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original  unsplit  panel 
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x— location  of  field  point 


Figure  4-4:  Investigation  of  the  error  due  to  representing  two  interior  source  panels 
with  a  single  panel.  The  length  of  the  unsplit  panel  is  1.0.  The  panel  is  split  in 
half  to  form  two  sub- panels.  The  panel  configuration  and  a  sketch  of  the  individual 
influences  are  shown  above  a  plot  of  the  ratio  •§£. 


weakly  on  the  initial  guess,  and  this  is  shown  in  Appendix  B. 

4.1.2  System  of  Equations 

For  a  given  cavity  planform,  Green’s  formula  should  be  satisfied  at  the  control 

points  of  all  N  x  M  —  ntPm  panels  on  the  key  blade.  As  before,  N  is  the  number 
of  panels  distributed  around  the  chord  of  the  blade  and  M  is  the  number  of  panels 
distributed  spanwise  (i.e.,  the  number  of  strips),  is  the  number  of  panels  on 
the  blade  which  are  “split”  by  the  trailing  edge  of  the  cavity  (Green’s  formula  is  not 
applied  on  the  split  panels).  In  addition  to  the  applications  of  Green’s  formula  on 
the  key  blade,  Green’s  formula  in  the  wake  is  satisfied  at  the  New  cavitating  panels 
on  Sew  ( New  does  not  count  the  split  panels  on  Sew )• 

The  final  system  of  equations,  including  the  influence  of  the  split  panels,  is  as 
follows: 


On  Scb 


M  ( N 
)  1  \  V,  Ainm<f>nm  4"  )  1  Ajnm^nm 

m— 1  l,  n=l  n=Nt,B+NcAVm+l 


Nia+NcAVm  9<j>  d(f> 

n=JVL,+ 1  01lnm  OTlL 


~  'y  ^  CilmQlm  ~~  Q 4"  W*lm  [(1  —  SOPm)^W»i 


-  <£lm]  |  = 


M  (  NcAVm  Nlb 

^  1  |  Vn+n)'T'^m  ’  ^nm  — 

msl  (  n=l  nsl 


y  1  Binm  &  ia«m  *  ~ 

n=^t*+^CXVm+l 


—  BimntrmU inmn,^  *  -^Ath  —  SOPmtVilm^ A’c>lVa,m  “  f  4 


+RHS5c 


t  =  1, ...,  N  x  M  —  n^. 


(4.10) 
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where 


NLB+NcAVm 

£im  =  ^  Ainm  +  Ai  +  SOPmWilm. 

n=iVi,a+ 1 


On  Sew 


M  (  Nib  N 
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—  53  FumQlm  -  QlmFimn.fnFLn  +  W^ilrn  [(1  —  SOPm)<^>iVm  —  ^lm]  /  — 
1=1  ) 


M  (  Ncav„ 

-4ir$i  +  <  -  5Z 

m=l  [  n=l 
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+  rhs5c 


*  =  1,...,JV  x  M 


(4.11) 


where 


WiJf+ATCilVm 

53  Anm  +  SOPmWilm  +  Din,^mFL„  +4 IT  if  *  is  on  the  mth  strip 

n=J^u+l 

NLB+NcAVm 

53  Anm  +  SOPmWilm  +  Din,^mFLm  otherwise 

n=Wx,a+l 


n«Pm  i»  the  chord  wise  index  of  the  split  panel  on  the  mth  strip. 

The  system  of  equations  (4.10)  and  (4.11)  must  be  solved  at  each  time  step.  This 
is  accomplished  by  an  accelerated  block-iterative  matrix  solver  [7,  17]. 


4.1.3  Cavity  Height  Computation 

Once  the  problem  has  been  solved  for  a  guessed  cavity  planform,  and  on  the 
cavity  panels  is  known,  the  cavity  height  (h(s,v,t)  on  the  blade  and  hw(s,  v,  t)  in  the 
wake)  can  be  determined  by  integrating  the  partial  differential  equations  (2.23)  and 
(2.30).  This  is  accomplished  by  replacing  the  partial  derivatives  of  h  and  hw  with 
two-point  backwards  difference  formulae  and  solving  for  h  and  hw  recursively.  Note 
that  the  derivatives  are  defined  at  the  control  points.  The  finite  difference  models  of 
the  derivatives  are  as  follows  ( hw  may  be  substituted  for  h): 


8h 

da 


ASi 


dh  ~  3 Jijm  - 

dv  3Avi  —  Av2 

dh  „  hL  -  /fc1 
dt  ~  At 


where 


hm  =  kh^  +  hl) 


Refer  to  Figure  4-5  for  clarification.  Substitution  of  the  finite  differences  in  (2.23)  and 
(2.30)  yields  recursive  formulae  for  hf+1  and  h*.+i  in  terms  of  previously  computed 
quantities.  The  height  of  the  cavity  at  its  trailing  edge,  S(r,t),  will  in  general  be 
non-zero,  unless  we  have  guessed  the  correct  cavity  planform.  The  means  by  which 
we  arrive  at  the  correct  planform  will  be  discussed  next. 


4.2  Finding  the  Cavity  Planform 

To  find  the  correct  discrete  cavity  planform  we  need  to  solve,  with  respect  to  the  M 
unknowns  1 1,  I3 f,  the  system  of  equations 

=  0  m  =  M  (4.12) 
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Figure  4-5:  Schematic  of  < 
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where  Sm  is  the  openness  of  the  cavity  trailing  edge  at  the  mth  spanwise  strip  and 
lm  is  the  cavity  length  at  midspan  of  the  same  strip.  The  system  of  equations  (4.12) 
may  be  solved  iteratively  by  applying  an  M-dimensional  Newton-Raphson  (secant 
method)  scheme.  The  updated  cavity  lengths  for  the  (k  +  l)rt  iteration  are  given  as 
follows: 


X*+1  =  -  J~'6k  +  Lk  (4.13) 

where 

LT  =  [hlih  ...  lm\ 

6T  =  [SiS,S3  ...  6m] 

and  where  J  is  the  Jacobian,  which  has  elements 

dSi 

*;  =  ~,  l<i,j<M. 

Ideally,  we  would  like  to  approximate  each  element  of  the  Jacobian  by  using  a 
two  point  finite  difference  scheme.  However,  this  would  mean  that  the  system  of 
equations  given  by  (4.10)  and  (4.11)  would  have  to  be  solved  M+l  times  for  each 
Newton-Raphson  iteration.  Even  in  the  case  of  an  accurate  initial  guess,  this  scheme 
clearly  would  be  prohibitively  expensive.  Fortunately,  the  iterative  method  may  be 
accelerated  by  applying  a  locally  scalar  Newton-Raphson  (secant  method)  scheme 
at  each  spanwise  strip,  which  is  equivalent  to  ignoring  the  off-diagonal  terms  of  the 
Jacobian  in  solving  (4.13)  for  the  updated  cavity  planform.  This  method  requires 
only  a  single  solution  of  the  system  of  equations  for  each  Newton-Raphson  iteration. 

The  second  guess  of  the  planform  is  taken  to  be  a  perturbation  of  the  initial 
guess  by  a  fixed  amount  (equal  to  2%  of  the  local  chord  length,  Cm)  in  a  direction 
determined  by  the  sign  of  the  local  6m  which  resulted  from  the  initial  guess.  For 
example,  if  Sm  from  the  first  guess  is  negative,  then  the  local  cavity  length  lm  is 
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shortened  by  .02cm  for  the  second  iteration.  In  general,  this  method  will  not  be 
stable  due  to  the  importance  of  the  off-diagonal  terms  of  the  Jacobian,  which  have 
been  ignored.  However,  numerical  damping  renders  the  method  stable  by  prohibiting 
each  cavity  length  lm  from  changing  by  an  amount  greater  than  a  multiple  of  Sm, 
|A/m|  <  |fe  *  5m|.  In  the  event  that  |£m|  is  large,  the  change  in  length  is  bounded  by 
|A/m|  <  0.05 Cm.  The  convergence  history  is  found  to  be  weakly  dependent  on  k,  and 
k  =  4  yields  satisfactory  convergence.  As  in  any  iterative  scheme,  however,  the  speed 
of  convergence  is  heavily  dependent  on  the  initial  guess.  Fortunately,  in  the  unsteady 
propeller  problem,  the  planform  from  the  previous  time  step  is  an  excellent  initial 
guess. 

At  each  time  step,  the  planform  will  be  considered  to  be  converged  when  the 
maximum  of  the  absolute  values  |5i| . . .  |5jvf|  is  less  than  or  equal  to  Suj.  The  sensitivity 
of  the  solution  to  6toi  is  discussed  in  section  4.2.2. 

To  demonstrate  the  three  dimensional  behavior  of  S  in  the  vicinity  of  the  correct 
cavity  planform,  Figure  4-6  shows  the  spanwise  distribution  of  5(y)  for  various  cavity 
planforms  on  a  rectangular  hydrofoil.  Among  these  cavity  planforms,  some  lie  to 
the  left  of  the  zero-5  planform,  and  some  lie  to  the  right.  Notice  that  5(y)  for  each 
planform  that  lies  to  the  left  of  the  correct  planform  is  everywhere  positive,  while  for 
those  that  lie  to  the  right  S(y)  is  everywhere  negative. 

Several  aspects  of  the  solution  are  best  discussed  before  describing  the  details  of 
the  unsteady  solution.  This  is  done  in  the  following  sections. 

4.2.1  Multiplicity  of  Solutions 

Figure  4-7  shows  a  series  of  planforms  for  a  rectangular  hydrofoil  at  an  angle  of  attack 
of  5°,  and  for  which  the  cavitation  number  is  varied  from  c  —  0.2  to  <r  =  1.6.  Below 
the  planforms,  the  nondimensional  cavity  length  at  the  midspan,  £(y  =  0)  is  plotted 
versus  the  ratio  f.  All  cavity  planforms  were  computed  using  the  same  discretization 
( N  =  80,  M  =  20).  In  order  to  avoid  clutter,  only  a  subset  of  the  cavity  planforms 
which  were  computed  and  presented  as  points  on  the  £(0)  vs  &  curve  have  been  shown 
in  the  upper  figure. 
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Figure  4-6:  Span  wise  £(y)  distributions  for  various  cavity  planforms  on  a  rectangular 
hydrofoil  which  has  a  NACA65a  crossection.  The  curves  are  labeled  1-10  to  correlate 
each  cavity  planform  to  its  6(y)  curve.  The  thickness  is  *  =  0.05  at  the  midspan, 
tapering  elliptically  to  zero  thickness  at  the  tip.  The  angle  of  attack  is  a  —  4°  and 
the  cavitation  number  is  a  =  0.4.  The  foil  is  twisted  spanwise,  and  the  *ng1»  of  twist 
is  given  by  at  =  -8y3  +  12yJ.  The  angle  of  twist  is  introduced  to  prevent  the  foil 
from  cavitating  near  the  tip. 
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The  most  striking  result  shown  in  Figure  4-7  is  the  apparent  multiplicity  of  cav¬ 
ity  planforms  for  a  given  cavitation  number.  Curves  #3  and  #4  both  correspond 
to  <t  =  0.9,  having  been  found  by  different  initial  guesses  (which  were  uniform  cav¬ 
ity  planforms  with  l(y)  =  0.9  and  1.3,  respectively).  Both  solutions  are  converged 
solutions  ( 6m  <  6tei  for  all  m). 

Except  in  a  certain  well  defined  regime,  there  is  a  one-to-one  correlation  between 
cavity  length  and  cavitation  number  for  two-dimensional  cavity  flows.  Inside  this 
regime,  the  multiplicity  of  solutions2  causes  no  difficulty,  since  we  know  which  solution 
is  the  physically  acceptable  one.  Although  it  has  generally  been  believed  that  the 
multiplicity  of  2-D  solutions  does  not  occur  for  3-D  solutions,  the  present  results 
appear  to  confirm  the  opposite.  In  fact,  the  £(0)  vs  ^  curve  looks  very  similar  to  the 
corresponding  2-D  result  [1,  66]. 

Note  that  each  solution  in  the  present  method  corresponds  to  a  “zero- crossing” 
of  S(y)  as  a  function  of  the  local  cavity  length  l(y).  In  our  solution,  only  the  zero- 
crossings  for  which  the  slopes  of  the  6(y)  vs  l(y)  curves  were  negative  were  searched 
for.  In  close  correspondence  with  the  multiplicity  which  occurs  in  two-dimensions,  our 
method  found  two  “negative-slope”  solutions.  Additional  solutions  may  exist  which 
our  numerical  scheme  does  not  search  for,  corresponding  to  §f(y)  >  0  (“positive- 
slope”)  or  to  “mixed-slope”  solutions  where  the  sign  of  ff  (y)  varies  with  y.  In  partic¬ 
ular,  for  the  high  aspect-ratio  rectangular  hydrofoil,  there  may  exist  a  “positive-slope” 
solution  which  lies  between  the  two  “negative-slope”  solutions.  However,  turning  off 
our  negative-slope  requirement  did  not  yield  a  positive-slope  solution  for  this  case. 
This  may  have  been  due  to  the  proximity  of  the  two  negative-slope  solutions.  Con¬ 
cerning  the  possible  existence  of  “mixed  slope”  solutions,  we  suspect  that,  if  they 
exist,  they  correspond  to  planforms  which  have  sharp  spanwise  discontinuities.  Our 
method  has  never  captured  such  a  solution.  It  is  clear  that  further  investigation, 
perhaps  including  a  cavitating  3-D  hydrofoil  experiment  to  examine  the  stability  of 
long  partial  cavities  and  short  supercavities,  is  necessary  to  resolve  this  issue. 

*For  some  values  of  the  parameter  *•  —  a  is  the  angle  of  attack  and  o  is  the  cavitation  n amber 
—  there  correspond  three  cavity  lengths,  bnt  only  the  shortest  is  physically  observed  [1]. 
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Figure  4-7:  Dependence  of  the  cavity  planform  on  <r,  showing  the  existence  of  multiple 
solutions. 
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As  would  be  expected,  the  multiplicity  vanishes  as  the  aspect  ratio  becomes 
smaller.  Figure  4-8  confirms  that  this  is  true  for  a  rectangular  foil  of  varying  span- 
to- chord  ratio. 


4.2.2  Sensitivity  of  the  Solution  to  6  Tolerance 

Figure  4-9  shows  the  dependence  of  the  computed  cavity  planforms  on  the  tolerance 
on  6.  The  planform  labeled  number  1  (which  is  also  the  converged  planform  for 
Std  =  0.1)  has  been  used  for  the  initial  guess  for  all  the  planforms  shown.  The  results 
shown  in  this  figure  suggest  that  Stoi  should  be  no  larger  than  .001. 

4.2.3  Nature  of  the  Tip  Solution 

Of  particular  importance  to  the  solution  is  the  behavior  of  the  tip  flow.  In  reality,  a 
cavitating  foil  or  propeller  blade  which  shows  significant  sheet  cavitation  typically  has 
a  cavitating  tip  vortex.  I  have  found  that  the  inviscid  solution  for  elliptic  hydrofoils, 
in  cases  where  the  tip  is  cavitating,  shows  a  tendency  toward  very  long  supercavities 
at  the  tip.  An  example  of  this  is  shown  in  Figure  4-10.  A  numerical  validation  for 
this  case  is  also  presented  in  Figure  5-10  in  section  5.2.2.  Due  to  high  skewness  of  the 
panels  near  the  tip,  it  is  difficult  to  obtain  accurate  representations  of  the  geometry  as 
well  as  accurate  influence  coefficients.  As  a  result,  it  is  often  difficult  for  the  solution 
to  converge  to  within  the  tolerance  set  for  6.  In  order  to  stop  the  iterations,  it  was 
sometimes  necessary  to  use  a  larger  tolerance  for  6  for  the  outer  1-2%  of  the  span. 
The  corresponding  6(y)  is  also  shown  in  Figure  4-10,  where  it  is  seen  that  S  is  largest 
at  the  tip. 

It  may  be  possible  to  improve  the  tip  solution,  perhaps  with  the  use  of  an  or¬ 
thogonal  grid  or  a  grid  which  conforms  to  the  local  streamlines,  but  this  was  not 
attempted  in  this  work.  One  reason  for  avoiding  these  refinements  is  the  apparent 
need  to  include  the  cavitating  tip  vortex  in  the  model.  It  will  be  proposed  that  the 
inner  sheet  cavity  solution  be  matched  to  an  outer  tip  vortex  solution,  in  which  the 
vortex  is  assumed  to  be  cavitating.  This  will  be  discussed  further  in  Chapter  6. 
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Figure  4-8:  Effect  of  aspect  ratio  on  multiplicity  of  solutions  for  a  rectangular  hydrofoil 
(i  =  .05, a  =  5°).  The  cavity  length  at  the  midspan  is  shown  vs  •  for  three  aspect 
ratios  AR  =  2, 5, 10. 
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Figure  4-9:  Dependence  of  the  cavity  pl&nform  on  6t*i,  for  a  rectangular  hydrofoil 
with  a  span-chord  ratio  of  5  at  an  angle  of  attack  of  5°.  (For  clarity,  the  span  and 
chord  are  not  drawn  to  scale.) 
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Figure  4-11:  Convergence  of  the  cavity  planform  with  updating  of  the  crossflow  terms. 
Same  foil  and  conditions  as  the  previous  figure. 

4.2.4  Convergence  of  the  Crossflow  Term 

It  was  mentioned  in  section  2.3  that  the  crossflow  terms  §£  and  V+  in  the  dynamic 
boundary  conditions  (2.17)  and  (2.20)  were  found  to  have  only  a  small  effect  on 
the  solution.  This  will  be  shown  for  three  geometries  in  this  section.  First,  Figure 
4-11  shows  the  cavity  planforms  from  two  consecutive  iterations  where,  in  between 
iterations  for  the  cavity  planform,  the  velocities  §£  and  V+  are  updated.  The  first 
term  is  computed  by  numerically  differentiating  the  potentials  using  second-order 
accurate  central  differences  (except  for  m  =  1  and  m  =  M,  where  forward  and 
backward  differences  are  used,  respectively).  The  second  term  is  computed  with 
central  differences  far  downstream  in  the  wake.  Since  the  wake  is  assumed  to  be 
force-free,  the  crossflow  velocity  is  constant  in  the  streamwise  direction  for  steady 
flow.  This  result  is  for  an  elliptic  foil  (the  same  foil  and  conditions  shown  in  Figure 
4-10).  Note  that  the  cavity  planforms  do  not  change  significantly  between  the  two 
iterations.  No  change  was  found  when  an  additional  iteration  was  tried.  From  this 
example,  which  showed  only  supercavitation,  it  may  be  concluded  only  that  the  term 
V+  has  little  effect  on  the  solution. 

To  gauge  the  importance  of  on  the  solution,  a  rectangular  hydrofoil  at  a  =  3° 
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Figure  4-12:  Convergence  of  the  cavity  planform  with  updating  of  the  crossflow  terms. 
Rectangular  hydrofoil  at  a  =  3°,  <r  =  0.5,  =  .014. 

is  tested  for  a  =  0.5  with  and  without  inclusion  of  the  crossflow  terms.  The  results 
are  shown  in  Figure  4-12.  Since  this  example  shows  only  partial  cavitation,  it  can  be 
seen  that  ||  also  has  a  negligible  effect  on  the  solution,  since  the  two  cavity  planforms 
are  nearly  indistinguishable. 

The  third  geometry  is  a  is  a  modified  AO- 177  propeller  (identical  to  the  AO- 177, 
except  that  the  skew  is  zero),  and  its  geometry  is  given  in  Table  4.1.  For  the  propeller 
solution,  the  crossflow  term  has  a  more  marked  effect  on  the  cavity  planform,  although 
the  convergence  is  still  very  quick.  Shown  in  Figure  4-13  are  cavity  planforms  on  the 
test  propeller  geometry  for  consecutive  crossflow  iterations. 

4.2.5  Nonlinear  Effect  of  Foil  Thickness 

The  nonlinear  effect  of  foil  thickness  on  cavity  extent  in  two-dimensional  partially 
cavitating  flows  is  well  known.  Namely,  increasing  the  foil  thickness  for  given  angle 
of  attack  and  cavitation  number  has  the  effect  of  decreasing  the  cavity  extent,  as  has 
been  shown  by  several  researchers  [34,  69,  72].  The  present  method  has  also  been 
found  to  manifest  the  nonlinear  thickness  effect  on  the  shape  of  the  cavity  for  2-D 
flows  [36].  A  similar  behavior  has  also  been  found  for  3-D  flows.  This  is  shown  in 


91 


r 

x 

— P — 
V 

rake 

TT _ 

skew(°) 

% 

hr 

0.20 

1.125 

Mil'll 

0.2070 

0.0490 

0.0414 

0.30 

1.223 

0.0058 

0.0 

0.2456 

0.0444 

0.0399 

0.40 

1.288 

0.0210 

0.0 

0.2722 

0.0367 

0.0361 

0.50 

1.318 

0.0410 

0.0 

0.2817 

0.0314 

0.0304 

0.60 

1.309 

0.0650 

0.0 

0.2684 

0.0300 

0.0236 

0.70 

1.250 

0.0913 

0.0 

0.2320 

0.0295 

0.0166 

0.80 

1.140 

0.1090 

0.0 

0.1815 

0.0281 

0.0107 

0.90 

0.970 

0.1168 

0.0 

0.1180 

0.0263 

0.0059 

0.95 

0.857 

0.1170 

0.0 

0.0810 

0.0252 

0.0038 

1.00 

0.722 

0.1148 

0.0 

0.0010 

0.0000 

0.0015 

Table  4.1:  Test  propeller  geometry,  identical  to  the  AO-177  propeller,  except  the  skew 
is  zero. 


— -  CROSSFLOW  NOT  INCLUDED 
-  CROSSFLOW  INCLUDED 


Figure  4-13:  Convergence  of  the  cavity  length  with  updating  of  the  crossflow  terms 
for  a  test  propeller. 


Figure  4-14:  Nonlinear  effect  of  foil  thickness  on  cavity  extent  and  volume.  Cavity 
pl&nforms  and  shapes  at  midspan  are  shown  for  rectangular  hydrofoils,  a  =  5°,  Aspect 
Ratio  AR  —  5.0,  t/c  =  .05,  .09  (span  and  chord  are  not  drawn  to  scale). 

Figure  4-14  where  the  cavity  planforms  and  shapes  at  midspan  are  shown  for  two 
rectangular  hydrofoils  at  two  cavitation  numbers.  Both  hydrofoils  have  an  aspect- 
ratio  of  5,  zero  camber,  and  a  NACA65a  thickness  section  the  maximum  thickness 
of  which  varies  elliptically  to  zero  in  the  spanwise  direction.  Both  foils  are  at  an 
angle  of  attack  a  =  5°.  The  first  has  a  thickness-to-chord  ratio  at  midspan  t/c  =  .05 
and  the  second  has  f/c  =  .09.  A  smaller  cavity  (in  volume  and  extent  at  all  spanwise 
locations)  is  predicted  for  the  thicker  foil  at  the  same  cavitation  number.  The  effect  of 
foil  thickness  on  cavity  size  appears  to  be  less  pronounced  in  the  case  of  supercavities. 
This  has  also  been  found  to  be  the  case  for  two  dimensional  flows. 
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4.3  The  Effect  of  0 o  on  the  Solution 

The  quantity  4>o  is  the  value  of  the  potential  at  the  leading  edge  of  the  cavity.  It  is 
not  a  primary  variable  because  it  is  defined  at  a  panel  edge  rather  than  at  the  control 
point.  Nonetheless,  it  is  an  unknown  in  the  system  of  equations  and  it  is  important 
to  the  solution.  Mainly,  it  is  important  for  <j> o  to  be  relatively  independent  of  gridding 
so  that  the  solution  will  converge  with  number  of  panels.  The  consequence  of  having 
a  grid  dependent  definition  is  shown  in  this  section. 

<f> o  was  first  defined  as  a  quadratic  extrapolation  of  the  unknown  potentials  at  the 
three  control  points  in  front  of  the  cavity  leading  edge  (see  Figure  3-5  for  clarifica¬ 
tion).  Figure  4-15  shows  the  resulting  convergence  characteristics  of  6(y)  for  a  3-D 
rectangular  hydrofoil  and  a  guess  of  the  cavity  planform.  S(y)  is  not  zero  because 
the  planform  is  not  correct.  The  rectangular  hydrofoil  is  5%  thick  at  the  midspan, 
tapering  elliptically  to  zero  thickness  at  the  tip.  The  angle  of  attack  is  a  =  3°  and  the 
cavity  length  is  £ (y)  =  0.5  for  all  y.  Computations  were  made  for  panel  arrangements 
of  (N,M)  =  (60, 30),  (80, 30),  (100, 30)  and  (120,30). 

The  convergence  of  S(y)  with  number  of  panels  is  seen  to  be  rather  slow.  This 
was  found  to  improve  drastically  when  a  cubic  extrapolation  is  used  (see  equation 
(3.11))  which  includes  a  condition  on  the  slope  of  <f>(s)  at  the  cavity  leading  edge. 
The  resulting  6(y)  convergence  is  shown  in  Figure  4-16. 


4.4  Time-Marching  Scheme 

The  time  marching  scheme  is  identical  to  that  used  in  the  fully  wetted  solution,  and 
is  described  in  detail  in  [22].  The  main  features  of  the  scheme  will  be  outlined  here 
for  the  sake  of  completeness. 

Time  is  discretized  into  equal  increments  of  span  At.  During  one  time  step,  each 
propeller  blade  rotates  through  an  incremental  angle  A6  =  wAt.  At  each  time  step, 
the  solution  is  found  for  the  key  blade  only,  while  the  singularities  on  the  other  blades 
are  assumed  to  be  known.  Before  proceeding  to  the  next  time  step,  vorticity  is  shed 
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Figure  4-15:  6(y)  for  a  rectangular  foil  at  a  =  4°  and  a  —  .16  with  a  cavity  length 
t(y)  =  0.5  for  all  y.  In  this  case  0o  is  defined  as  a  quadratic  extrapolation. 


Figure  4-16:  S(y)  for  a  rectangular  foil  at  a  =  4°  and  a  =  .16  with  a  cavity  length 
l(y)  =  0.5  for  all  y.  In  this  case,  <f> 0  is  defined  as  a  cubic  extrapolation  with  specified 
slope  at  the  cavity  leading  edge. 


downstream  along  the  assumed  wake  surface  through  an  angular  distance  equal  to 
the  incremental  rotation  A0.  This  enforces  the  equality  of  the  strength  of  the  shed 
vorticity  and  the  time  rate  of  change  of  the  circulation  on  the  blade.  The  strength 
of  the  circulation  is  given  in  terms  of  the  potential  jump  at  the  trailing  edge  of  the 
blade 

r(r,t)  =  A  =  <f>j  (r,t)  =  <f>sm 

and  the  potential  jump  at  the  first  wake  panel  is  given  by 


A  <j>(r,t) 


T(i)  +  T(t  -  At) 

2 


(4.14) 


Thus,  the  vorticity  convection  is  used  to  define  the  wake  dipole  strengths.  Although 
this  an  extension  of  Morino’s  Kutta  condition  in  steady  flow,  it  is  not  equivalent  to 
the  Kutta  condition  he  applied  in  unsteady  flow  [52]. 

The  solution  is  initiated  by  the  fully  wetted  steady  solution.  Next,  the  fully  wetted 
unsteady  solution  is  obtained,  which  then  serves  as  the  “initial  guess”  for  the  unsteady 
cavity  solution.  The  unsteady  cavity  solution  is  turned  on  when  the  key  blade  is  at 
the  6  O’Clock  position,  so  that  it  is  likely  to  start  out  fully  wetted.  This  was  found 
to  be  an  important  time-saving  measure  in  the  linear  solution  by  Lee  [45].  The  singu¬ 
larity  strengths  on  the  other  blades  and  their  wakes  are  taken  from  earlier  key  blade 
solutions.  During  the  first  revolution  of  the  key  blade,  the  singularities  on  another 
blade  are  taken  either  from  the  fully  wetted  solution  or  from  the  cavity  solution  when 
the  key  blade  was  in  the  same  angular  position.  For  subsequent  revolutions,  the  other 
blade  singularities  are  taken  from  previous  cavity  solutions  when  the  key  blade  was 
in  the  same  angular  position. 


4.5  Treatment  of  the  Unsteady  Terms 

The  dynamic  and  kinematic  boundary  conditions  (2.17),  (2.19),  (2.23)  and  (2.30),  in¬ 
clude  time  derivatives  of  the  potential  and  cavity  thickness.  The  numerical  treatment 
of  these  terms  will  be  described  in  this  section. 
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The  dynamic  boundary  condition  (2.17)  expresses  <f>  as  a  function  of,  among  other 
things,  the  time  derivative  §£.  At  a  given  time  step,  §&  is  assumed  to  be  known  and 
given  by  a  quadratic  extrapolation  of  its  values  from  previous  time  steps.  Specifically, 
the  derivative  is  assumed  to  behave  locally  quadratically  in  time 

^(t)  ~  At3  +  Bt  +  C 

at 

and  the  coefficients  A,BfC  are  defined  in  terms  of  the  values  of  at  the  three 
previous  time  steps.  The  time  derivative  at  the  current  time  step  is  then  extrapolated 
forward  according  to 

w<‘>  =  3^‘  -  A‘>  -  3w((  -  2A‘>  +  8W<‘  - 3Ai) 

The  derivative  at  each  time  step  is  computed  numerically  via  a  backwards  fourth- 
order  accurate  finite  difference  scheme3.  The  extrapolation  of  §|  fits  well  with  the 
time  marching  scheme,  since  several  revolutions  are  required  to  obtain  the  steady- 
state  oscillatory  solution  due  to  the  time  lag  of  the  solution  on  the  other  blades.  The 
convergence  of  the  solution  with  revolutions  is  discussed  in  section  5.1. 

In  the  kinematic  boundary  conditions  (2.23)  and  (2.30),  the  time  derivatives  of 
the  cavity  thickness,  and  are  replaced  by  first  order  backwards  finite  differ¬ 
ences  (see  section  4.1.3).  These  do  not  require  any  further  special  treatment.  These 
boundary  conditions  are  used  to  compute  the  cavity  thickness,  after  the  solution  at 
the  current  time  step  has  been  found  and  the  source  strengths  are  known.  In  addition, 
however,  these  terms  must  be  used  to  determine  the  strength  of  a  source  panel  which 
is  wetted  at  the  current  time  step,  but  which  was  cavitating  during  the  previous  time 
step.  In  other  words,  a  collapsing  cavity  has  an  additional  sink  on  the  panels  which 
lie  between  the  trailing  edge  of  the  cavities  from  consecutive  time  steps.  This  was 

’Various  finite  difference  schemes  have  been  tried  for  the  computation  of  the  fully  wetted  unsteady 
pressure  distribution  and  the  fourth  order  finite  difference  model  was  chosen  for  reasons  of  robustness 
[22]. 
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first  noticed  by  Van  Houten  [73].  The  magnitude  of  the  source  will  be  given  as 


8£ 

dn 


(4.15) 


In  order  to  understand  the  importance  of  the  unsteady  terms  of  the  boundary 
conditions,  an  unsteady  two  dimensional  model  of  a  partially  cavitating  hydrofoil  was 
investigated,  wherein  the  cavitation  number  is  allowed  to  vary  harmonically  (<r(t)  = 
<r0  +  <rg  cosu ft)  while  the  inflow  is  steady  and  uniform  and  the  wake  is  steady  (no  shed 
vorticity,  but  the  wake  dipole  strength  changes  each  time  step).  Although  this  model 
is  not  an  adequate  model  of  the  physical  flow,  it  allows  us  to  see  the  relative  effects  of 
the  various  unsteady  terms  on  the  global  2-D  solution.  To  this  end,  the  time  history  of 
the  cavity  length  was  computed  for  a  NACA16006  foil  operating  at  an  angle  of  attack 
of  4°  and  with  a  time-dependent  cavitation  number  equal  to  <r(t)  =  1.2  +  0.2  cos  art 
at  a  reduced  frequency  of  k  =  =  1.0.  Figure  4-17  shows  the  cavity  length  as  a 

function  of  u>t  for  one  cycle  for  two  different  solutions.  The  first  solution  contains  no 
history  at  all;  the  time  derivatives  of  both  the  potential  and  the  cavity  thickness  are 
set  to  zero.  The  second  solution  contains  the  unsteady  terms.  Note  that  the  unsteady 
solution  is  shifted  to  the  right,  indicating  that  the  growth  stage  lasts  longer  than  the 
collapse  stage,  a  well  known  characteristic  of  unsteady  cavitation  [73,  66]. 

Figure  4-18  shows  cavity  shapes  from  several  time  steps  during  the  second  cycle 
of  the  unsteady  solution.  Of  particular  interest  is  the  difference  between  the  cavity 
shapes  at  u>t  =  90°  and  art  =  270°  for  which  the  cavitation  numbert  are  identical  The 
difference  between  the  shapes  is  due  entirely  to  the  unsteady  terms.  This  shows  the 
importance  of  the  unsteady  terms. 
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Figure  4-17:  Cavity  length  vs  uit  for  a  NACA16006  2-D  section  at  a  =  4°  and 
<r(t)  =  1.2  +  0.2 cosuit.  Shown  are  the  quasi-steady  and  the  unsteady  solutions. 


Figure  4-18:  Cavity  shapes  from  several  time  steps  of  the  unsteady  solution  shown  in 
the  previous  figure. 
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Chapter  5 


Numerical  Validation 


The  goal  of  this  chapter  is  to  show  the  present  method  to  be  convergent  and  consistent, 
as  all  viable  numerical  methods  must  be.  Here,  a  convergent  solution  implies  that  as 
the  number  of  panels  tends  to  infinity  and  the  average  panel  size  tends  to  zero  the 
solution  would  converge  to  an  answer  which  is  exact  for  the  given  assumptions.  A 
consistency  test  is  applied  to  test  the  accuracy  of  the  solution.  Comparison  is  made 
between  the  results  of  the  present  theory  with  results  of  linear  theory  and  linear  theory 
with  leading  edge  corrections.  Some  comparisons  are  made  to  published  experimental 
results  and  observations. 


5 . 1  Convergence 

Figure  5-1  shows  convergence  of  the  spanwise  S  distribution  and  chordwise  cavity 
thicknesses  for  a  rectangular  hydrofoil  and  a  guess  of  the  cavity  planform  l  =  0.5  at 
all  spanwise  strips.  This  was  the  initial  convergence  test  and  eliminates  the  effects  of 
the  split-panel  method  ( l  =  0.5  always  corresponds  to  a  panel  boundary  with  cosine 
spaced  panels)  and  the  5-tolerance  associated  with  the  planform  convergence  scheme. 
The  hydrofoil  has  a  NACA65a  crossection  with  *  =  .05  at  the  midspan,  tapering 
elliptically  to  zero  at  the  foil  tips.  The  foil  is  twisted  so  that  the  angle  of  attack 
is  lower  at  the  tips  than  at  the  midspan;  the  twist  angle,  in  degrees,  is  given  by 
at(y)  =  8y3  -  12y*.  The  detachment  point  is  2.4%  of  the  local  chord  length  aft  of 
the  leading  edge  at  each  spanwise  strip.  The  aspect  ratio  is  AR  =  5.9,  the  cavitation 
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number  is  <x  =  0.4  and  the  angle  of  attack  is  a  =  4°.  Figure  5-1  displays  excellent 
convergence  characteristics  for  these  variables  under  these  conditions. 

Figure  5-2  shows  the  convergence  of  the  cavity  planform  and  the  cavity  volume 
with  number  of  panels  for  the  same  rectangular  hydrofoil  used  in  the  previous  exam¬ 
ple,  except  that  in  this  case  the  twist  angle  at  is  zero  everywhere.  This  convergence 
test  is  more  demanding  than  the  previous  because  of  the  use  of  the  split  panel  method, 
which  introduces  some  error  into  the  computation.  The  angle  of  attack  is  3°  and  the 
cavitation  number  is  0.12.  The  chordwise  discretization  is  varied  from  60  to  120  with 
the  spanwise  discretization  fixed  at  20.  With  the  chordwise  discretization  fixed  at  80, 
the  spanwise  discretization  is  varied  from  20  to  30.  In  each  case,  the  wake  is  divided 
into  60  panels,  with  40  in  the  near  wake,  which  extends  one  half-span  downstream, 
and  20  in  the  far  wake,  which  extends  an  additional  9  half-spans.  From  these  re¬ 
sults,  it  can  be  concluded  that  a  foil  discretization  of  80  x  20  is  sufficient  to  obtain  a 
converged  result  for  the  rectangular  hydrofoil. 

The  next  task  is  to  test  the  convergence  characteristics  of  the  propeller  solution. 
Figure  5-3  shows  convergence  of  the  cavity  planform  and  cavity  thickness  on  a  test 
propeller  which  is  a  N4381  geometry  with  modified  pitch  distribution.  The  modified 
N4381  has  an  unloaded  tip  in  order  to  ensure  that  that  the  blades  do  not  supercavitate 
near  the  tip,  so  that  the  convergence  characteristics  of  the  partial  cavitation  solution 
may  be  investigated.  The  propeller  geometry  is  given  in  Table  5.1.  First,  the  method 
is  applied  to  the  test  propeller  with  a  single  blade  operating  in  a  uniform  inflow  with 
advance  coefficient  Js  =  Vs/(nD)  =  0.8,  where  Vs  is  the  ship  speed,  and  <rn  =  2.7. 
The  variation  in  hydrostatic  pressure  is  ignored  and  the  resulting  solution  is  steady. 
The  cavity  is  assumed  to  detach  2.4%  of  the  local  chord  length  aft  of  the  leading 
edge.  Results  of  this  convergence  test  show  that  a  blade  discretization  of  80  x  20  is 
sufficient  to  obtain  a  converged  result  in  steady  flow  for  the  partially  cavitating  test 
propeller. 

The  convergence  of  a  supercavitating  propeller  is  tested  next.  Shown  in  Figure 
5-4  is  the  steady  cavity  solution  on  a  single  bladed  version  of  the  AO- 177  propeller 
at  Js  =  0.6  and  «rn  =  2.5.  The  AO- 177  propeller  geometry  is  provided  in  Table  5.2. 
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Figure  5-1:  Convergence  of  span  wise  and  chord  wise  cavity  thickness  distributions 
with  number  of  panels.  The  foil  and  cavity  characteristics  are  described  in  the  text. 
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Figure  5-2:  Convergence  of  cavity  planform  (top  plot)  and  cavity  volume  (bottom 
plot)  with  number  of  panels. 
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Figure  5-3:  Convergence  of  cavity  planform  and  cavity  thicknesses  with  number  of 
panels  for  the  modified  N4381  one-bladed  propeller  in  steady  flow;  J$  =  0.8 ,  an  =  2.7. 

The  cavity  detaches  2.4%  of  the  local  chord  length  aft  of  the  leading  edge  at  all 
span  wise  locations. 
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F 

K 

* 

ia 

mm 

hr 

0.20 

1.332 

0.0 

0.0 

0.174 

0.0351 

0.0651 

0.25 

1.338 

0.0 

0.0 

0.202 

0.0369 

0.0594 

0.30 

1.345 

0.0 

0.0 

0.229 

0.0368 

0.0537 

0.40 

1.358 

0.0 

0.0 

0.275 

0.0348 

0.0441 

0.50 

1.350 

0.0 

0.0 

0.312 

0.0307 

0.0360 

0.60 

1.310 

0.0 

0.0 

0.337 

0.0245 

0.0286 

0.70 

1.250 

0.0 

0.0 

0.347 

0.0191 

0.0219 

0.80 

1.150 

0.0 

0.0 

0.334 

0.0148 

0.0157 

0.90 

0.950 

0.0 

0.0 

0.280 

0.0123 

0.0101 

0.95 

0.750 

0.0 

0.0 

0.210 

0.0128 

0.0072 

1.00 

0.500 

0.0 

0.0 

0.000 

0.0120 

0.0044 

Table  5.1:  The  DTMB  propeller  N4381,  with  modified  pitch  distribution  to  unload 
the  tip  (referred  to  in  the  text  as  the  “test  propeller”). 


In  this  case,  the  cavity  is  assumed  to  detach  at  the  blade  leading  edge.  To  find  the 
planform,  the  trailing  edge  openness  6  is  required  to  converge  to  zero  everywhere 
except  at  the  two  outermost  strips,  for  r/R  >  .98.  As  mentioned,  the  converged 
planform  is  difficult  to  obtain  near  the  tip  due  to  inaccurate  modelling  of  the  geometry 
(see  discussion  in  section  4.2.3).  However,  in  most  cases  S(r)  converges  to  less  than 
10%  of  the  tip  chord  length  (which  is  usually  less  than  5%  of  the  blade  diameter).  In 
contrast,  the  tolerance  set  for  Sm  for  r/R  <  .98  in  this  example  was  6toi  =  -001. 

Next,  the  test  propeller  is  run  in  a  nonuniform  inflow  which  is  characterized  by 
a  15%  dent  in  the  axial  inflow,  symmetric  about  9  =  0  (see  equation  (2.2)  for  a 
definition  of  the  inflow  velocity).  The  other  components  of  the  flow  are  zero.  Due  to 
the  circumferential  symmetry,  all  of  the  harmonic  coefficients  of  the  sine  series  for  the 
axial  flow  are  zero  ( Bn  —  0  Vn).  The  mean  axial  velocity  is  set  to  one  (A\(r)  =  1). 
The  coefficients  of  the  cosine  series  are  given  in  Table  5.3. 

The  circumferential  distribution  of  axial  inflow  velocity  is  shown  in  Figure  5-5. 
Also  shown  in  this  figure  are  convergence  with  revolutions  and  convergence  with  time 
step  size  of  the  cavity  volume  history. 
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T 

: h 

p 

T5 

skew(°) 

C 

T5 

*** 

hr 

0.20 

1.125 

minima 

0.2070 

0.0490 

0.0414 

0.30 

1.223 

0.0058 

2.2 

0.2456 

0.0444 

0.0399 

0.40 

1.288 

0.0210 

7.1 

0.2722 

0.0367 

0.0361 

0.50 

1.318 

0.0410 

13.1 

0.2817 

0.0314 

0.0304 

0.60 

1.309 

0.0650 

20.0 

0.2684 

0.0300 

0.0236 

0.70 

1.250 

0.0913 

27.7 

0.2320 

0.0295 

0.0166 

0.80 

1.140 

0.1090 

34.5 

0.1815 

0.0281 

0.0107 

0.90 

0.970 

0.1168 

40.3 

0.1180 

0.0263 

0.0059 

0.95 

0.857 

0.1170 

42.79 

0.0810 

0.0252 

0.0038 

1.00 

0.722 

0.1148 

45.0 

0.0010 

0.0000 

0.0015 

Table  5.2:  The  AO- 177  propeller  geometry. 


n 

An 

n 

An 

1 

1.0 

9 

0.0000 

2 

-.0360 

10 

0.0009 

3 

-.0318 

11 

0.0009 

4 

-.0257 

12 

0.0005 

5 

-.0187 

13 

0.0000 

6 

-.0120 

14 

-.0003 

7 

-.0064 

15 

-.0003 

8 

-.0023 

16 

-.0002 

Table  5.3:  Harmonic  coefficients  An  for  the  nonuniform  axial  inflow  velocity. 


Figure  5-4:  Convergence  of  cavity  planform  with  number  of  panels  for  the  one-bladed 
AO-177  propeller  in  steady  flow;  Js  =  0.6,  <rn  =  2.5.  Shown  also  is  the  discretized 
blade  for  ( N,M )  =  (100,20). 
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Figure  5-5:  Convergence,  with  revolutions  and  with  time  step  size,  of  cavity  volume 
history  for  the  one-bladed  test  propeller  operating  in  a  nonuniform  axial  inflow.  The 
axial  inflow  is  shown  at  the  top  of  the  figure. 
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5.2 


Consistency 


5.2.1  Fully  Wetted  Consistency  Test 


One  very  useful  consistency  test  consists  of  replacing  the  dynamic  boundary  condition 
(3.9)  with 


4>(NLB+n)m  =  <t>fw 


n  =  1,  Ncxvm  +  Ncm 
m  =  1, M 


(5.1) 


where  <f>fw  is  the  potential  from  the  fully  wetted  solution.  The  resulting  solution 
should  be  equivalent  to  the  fully  wetted  solution  (i.e.,  the  cavity  thickness  should 
be  zero  everywhere  and  the  resulting  circulation  distribution  should  be  equivalent  to 
the  fully  wetted  circulation).  As  a  test  of  self-consistency,  this  tool  was  very  useful 
for  de-bugging  the  computer  code.  The  first  check  of  any  code  was  to  apply  this 
consistency  test. 

Figure  5-6  shows  the  results  of  the  fully  wetted  consistency  test  for  the  N4381  test 
propeller  in  steady  flow.  Shown  in  the  figure  are  the  circulation  distributions  from  the 
fully  wetted  solution  and  the  cavity  solution  (with  the  dynamic  boundary  condition 
replaced  with  (5.1)).  The  two  circulation  distributions  are  very  nearly  identical,  with 
the  slight  discrepancy  near  the  tip  perhaps  attributable  to  the  different  tolerances  set 
on  the  iterative  matrix  solutions.  This  discrepancy  is  considered  negligible. 


5.2.2  Pressure  Validation 

To  validate  the  solution,  a  modified  geometry  is  created  by  combining  the  cavity 
and  foil  (or  blade)  to  form  a  single  surface  and  representing  the  resulting  surface 
by  a  lattice  of  quadrilateral  panels.  The  grid  is  adapted  to  the  shape  of  the  cavity 
so  that  panels  are  clustered  around  the  cavity  trailing  edge.  The  fully  wetted  flow 
is  then  computed  about  the  modified  foil  and  the  resulting  pressure  distribution  is 
compared  to  the  pressure  from  the  cavity  solution.  If  the  cavity  shape  is  the  correct 
shape  which  corresponds  to  the  applied  dynamic  boundary  condition,  then  the  two 
pressure  distributions  will  match.  For  example,  if  the  solution  were  the  fully  nonlinear 
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Figure  5-6:  Results  of  a  fully  wetted  consistency  test.  Circulation  distributions  from 
the  fully  wetted  solution  and  the  cavity  solution.  The  two  distributions  should  match 
because  the  potential  on  the  cavity  panels  is  set  equal  to  the  fully  wetted  potential. 

solution,  then  the  pressures  should  match  exactly.  The  discrepancy  between  the  two 
pressure  distributions  is  therefore  a  measure  of  the  error  introduced  by  applying  the 
boundary  conditions  on  the  approximate  surface.  This  test  will  be  referred  to  as 
pressure  validation ,  and  will  be  shown  for  several  configurations. 

3-D  Partially  Cavitating  Hydrofoil 

Figures  5-8  and  5-9  show  the  results  of  a  pressure  validation  test  for  a  partially 
cavitating  rectangular  hydrofoil  which  is  the  same  as  the  foil  shown  in  Figure  5-1 
(including  the  angle  of  twist  which  tends  to  eliminate  tip  cavitation).  An  example  of 
an  adapted  grid  on  which  the  fully  wetted  flow  is  computed  is  shown  in  Figure  5-7. 
In  Figure  5-8  the  angle  of  attack  is  a  =  3°  and  ~  =  .105,  while  in  Figure  5-9  the  angle 
of  attack  is  a  =  5°  for  the  same  For  the  smaller  of  the  two  angles  of  attack,  the 
comparison  between  the  fully  wetted  pressure  and  the  cavity  pressure  is  excellent  at 
all  spanwise  strips.  For  the  higher  angle,  however,  the  comparison  deteriorates.  This 
indicates  that  the  shape  of  the  cavity  horn  the  first  iteration  is  further  from  the  correct 
nonlinear  shape  for  higher  angles  of  attack.  This  is  not  surprising.  However,  the 
results  in  2-D,  discussed  extensively  in  Appendix  B,  show  that  the  difference  between 
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the  first  iteration  and  converged  cavity  lengths  and  volumes  is  not  substantial,  even 
at  moderately  high  angles  of  attack.  In  other  words,  from  my  experience  in  applying 
the  same  validation  test  in  2-D,  the  pressure  distributions  from  the  first  iteration  and 
the  converged  solution  differ  more  than  the  corresponding  cavity  lengths  and  shapes. 


3-D  Supercavitating  Hydrofoil 

Figure  5-10  shows  pressure  validation  for  an  elliptic  hydrofoil  (the  chord  length  is 
maximum  at  the  midspan,  tapering  elliptically  to  zero  at  the  tip)  at  a  =  3°  and 
f  =  0.33.  The  aspect  ratio  is  AR  =  1.27.  In  general,  this  comparison  is  very  good 
everywhere  except  very  close  to  the  tip.  For  the  foil  strips  beyond  y  =  .98R,  the 
planform  was  not  required  to  converge  to  zero-6,  although  it  came  very  close.  The 
difficulty  in  finding  the  correct  planform  near  the  tip  would  be  mitigated  by  the 
inclusion  of  the  tip  vortex  in  the  model.  This  will  be  discussed  further  in  Chapter  6. 

One-bladed  Test  Propeller  in  Steady  Flow 

Figure  5-11  shows  pressure  validation  for  the  test  propeller  in  steady  flow1  at  r/R  = 
.75.  The  advance  coefficient  is  Js  =  0.8  and  the  cavitation  number  is  <r„  =  2.7. 
For  these  conditions,  the  pressure  validation  shows  that  the  first  iteration  propeller 
solution  does  a  good  job  of  predicting  the  correct  nonlinear  cavity  shape  at  this 
advance  coefficient. 

5.2.3  Comparison  to  PUF-3A 

Figure  5-12  shows  a  comparison  between  the  cavity  planform  computed  by  the  present 
method  (labeled  PROPCAV,  which  is  the  name  of  the  program)  and  those  computed 
by  linear  theory  and  linear  theory  with  leading  edge  corrections  (labeled  PUF-3A 
[31,  34]).  The  computations  are  done  for  steady  flow  conditions  on  the  one-bladed 
modified  N4381  test  propeller.  Linear  theory  is  seen  to  overpredict  the  cavity  extent, 

'The  pressure  validation  test  is  not  applied  to  the  unsteady  case  due  to  the  excessive  computation 
times  that  would  require. 
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Figure  5-7:  Adapted  grid  fit  to  the  shape  of  the  cavity,  used  to  compute  the  fully 
wetted  flow  on  the  modified  foil/cavity  boundary  for  the  pressure  validation  shown 
in  the  next  figure. 


Figure  5-8:  Pressure  validation  for  a  rectangular  hydrofoil  at  a  =  3°  and  *  =  .105. 
The  cavity  detaches  at  .4%  of  the  chord  length  aft  of  the  leading  edge  at  all  spanwise 
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Figure  5-9:  Pressure  validation  for  a  rectangular  hydrofoil  at  a  =  5°  and  ®  =  .105. 
The  cavity  detaches  at  .4%  of  the  chord  length  aft  of  the  leading  edge  at  all  spanwise 
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Figure  5-10:  Pressure  validation  for  a  supercavitating  elliptic  hydrofoil  at  a  =  3°  and 
«  =  .33. 
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Figure  5-11:  Pressure  validation  for  the  test  propeller  at  Js  =  0.8  and  a  =  2.7. 

as  it  does  in  2-D.  We  see  that  the  cavity  extent  is  overpredicted  by  linear  theory  in 
both  spanwise  and  chordwise  directions.  On  the  other  hand,  the  linear  theory  with 
leading  edge  corrections  is  seen  to  underpredict  the  extent,  although  the  shape  it 
produces  is  more  accurate  than  linear  theory.  A  comparison  is  then  made  for  the 
same  propeller  in  nonuniform  axial  wake  inflow  (the  same  inflow  as  shown  in  Figure 
5-5).  In  this  case,  the  hydrostatic  terms  are  turned  on  ( Fr  —  n2D/g  =  10.45).  The 
advance  coefficient  and  the  cavitation  number  are  kept  the  same  as  in  the  steady  flow 
case.  The  cavity  volume  histories  predicted  by  the  three  methods  are  shown  in  Figure 
5-13. 


5.2.4  Effect  of  the  Hub 

The  hub  may  be  included  in  the  solution.  Quadrilateral  panels  are  placed  on  the 
surface  of  the  hub  and  Green’s  formula  is  satisfied,  subject  to  the  kinematic  boundary 
condition.  The  modelling  of  the  hub  has  been  described  by  J.T.  Lee  [49]. 

The  effect  of  the  hub  on  the  cavity  solution  is  shown  in  Figure  5-14  for  the  test 
propeller  at  Js  =  0.8  and  <rn  =  2.7  in  steady  flow.  Note  that  the  presence  of  the  hub 
makes  the  cavity  larger  at  the  inner  radii  because  the  local  increase  in  circulation. 
This  effect  is  not  captured  in  the  linear  solution  (PUF-3A),  because  the  hub  is  not 
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Figure  5-12:  Cavity  planforms  predicted  by  PUF-3A  (with  and  without  the  leading 
edge  correction)  and  the  present  method  (implemented  in  the  code  PROPCAV)  on 
the  modified  N4381  test  propeller  at  Js  =  0.8  and  <r  =  2.7  in  uniform  flow. 
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Figure  5-13:  Cavity  volume  histories  predicted  by  PUF-3A  (with  and  without  the 
leading  edge  correction)  and  the  present  method  on  the  test  propeller  at  Js  —  0.8 
and  <r  =  2.7  in  nonuniform  flow. 
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included  in  that  model. 


5.2.5  Comparison  to  Experimental  Results 

No  systematic  comparisons  of  numerical  results  to  experimental  results  have  been 
done  for  this  thesis.  However,  several  qualitative  comparisons  can  be  made  to  pub¬ 
lished  observations.  Since  one  of  the  main  accomplishments  of  the  present  work  has 
been  the  prediction  of  cavity  planforms,  it  would  be  useful  to  compare  predicted 
planforms  to  observed  planforms.  In  a  paper  by  Kato,  et  al  [26],  experimental  re¬ 
sults  are  presented  for  an  elliptic  hydrofoil.  A  reproduction  of  one  of  their  observed 
planforms  is  shown  in  Figure  5-15.  Although  there  are  some  ambiguities  about  the 
geometry  and  operating  conditions,  I  have  computed  a  planform  using  the  present 
method  for  conditions  which  I  believe  are  close  (if  not  identical)  to  the  experimental 
conditions.  The  computed  results  are  also  shown  in  Figure  5-15.  This  comparison 
shows  at  least  qualitative  agreement.  For  example,  the  behavior  of  the  cavity  near 
the  tip  (for  r/R  >  .8,  say)  seems  to  agree  qualitatively. 

Van  Houten  studied  a  3-D  pitching  rectangular  foil  under  cavitating  and  noncavi- 
tating  conditions  [72].  The  unsteady  nature  of  his  experiment  makes  it  impossible  to 
compare  his  observations  to  the  present  computations.  However,  the  rapid  decrease 
of  the  cavity  extent  near  the  tip,  which  may  be  expected  to  occur  also  in  steady  flow, 
appear  in  the  present  computations  (see,  for  example,  Figure  5-2).  An  example  of 
Van  Houten’s  experimental  observations  is  shown  in  Figure  5-16,  where  this  sharp 
drop-off  of  cavity  length  near  the  tip  is  visible. 

Cumberbatch  studied  tip  effects  on  a  cavitating  rectangular  wedge  at  an  angle 
of  attack.  In  his  study,  he  observed  supercavity  patterns  which  tended  to  disappear 
near  the  tip,  in  qualitative  agreement  with  the  present  results  (as  shown  in  Figure 
5-2).  According  to  Cumberbatch’s  observations,  the  tip  vortex  and  the  sheet  cavity 
were  separated  by  a  wetted  region. 

A  supercavitating  2-D  hydrofoil  experiment  was  carried  out  at  MIT’s  Marine 
Hydrodynamics  Laboratory  in  the  Variable  Pressure  Water  Tunnel  by  Kinnas  and 
Mazel  [40].  In  that  work,  the  authors  compared  measured  velocities  to  those  computed 
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experimental 


Figure  5-16:  A  reproduction  of  a  photograph  of  a  partially  cavitating  3-D  pitching 
hydrofoil.  The  cavitation  number  was  a  —  .97  and  the  reduced  frequency  of  the 
pitching  motion  was  k  =  1.0.  The  picture  was  taken  by  Van  Houten  and  found  in 
MIT  Marine  Hydrodynamics  Lab  archives. 

by  the  fully  nonlinear  boundary  element  method  developed  in  this  thesis.  In  order  to 
make  the  comparison,  the  method  of  images  was  used  to  account  for  the  presence  of 
the  walls.  An  example  of  the  resulting  comparison  is  shown  in  Figure  5-17.  Shown 
are  several  numerical  and  experimental  horizontal  velocity  profiles  computed  and 
measured  near  the  surface  of  the  cavity.  The  comparison  is  very  good  away  from  the 
cavity  surface,  deteriorating  as  the  measurement  point  approaches  the  mean  cavity 
surface.  One  explanation  for  the  poor  comparison  near  the  cavity  surface  may  be 
the  transient  nature  of  the  flow  there  due  to  unsteady  motions  of  the  cavity  [40]. 
Nevertheless,  the  authors  concluded  that  the  overall  agreement  between  theory  and 
experiment  was  very  good,  with  the  analysis  predicting  the  flow  away  from  the  cavity 
quite  accurately. 
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Figure  5-18:  A  comparison  of  experimental  (measurements  by  Kinnas  and  Mazel)  and 
numerical  (computations  by  Fine)  velocity  profiles  on  a  supercavitating  2-D  hydrofoil 
in  a  water  tunnel.  Taken  from  Kinnas  and  Mazel  (1992). 


Chapter  6 


Conclusions  and  Recommendations 

6.1  Conclusions 

In  this  thesis,  an  efficient  and  robust  numerical  method  has  been  developed  for  accu¬ 
rately  predicting  the  extent  and  shape  of  sheet  cavities  on  propellers  in  unsteady  flow. 
Along  the  way,  the  method  was  extensively  validated  by  numerical  convergence  and 
consistency  tests  for  2-D  and  3-D  hydrofoils  in  steady  flow,  for  propellers  in  steady 
flow,  and  for  propellers  in  unsteady  flow.  The  applied  boundary  conditions  include 
all  nonlinear  and  crossflow  terms.  However,  for  three  dimensional  flows,  the  condi¬ 
tions  are  satisfied  on  the  solid  body  surface  (either  the  hydrofoil  or  propeller  blade 
surface)  and  the  portion  of  the  wake  surface  overlapped  by  the  cavity,  rather  than  on 
the  exact  cavity  surface.  This  solution  represents  the  first  iteration  towards  the  fully 
nonlinear  solution.  The  first  iteration  solution  has  been  shown  (here  and  in  the  cited 
references),  for  two-dimensional  flows,  to  be  very  close  to  the  converged  nonlinear 
solution  for  a  reasonable  range  of  angles  of  attack.  This  result  has  been  used  to  argue 
that  the  first  iteration  is  sufficient  to  obtain  a  nonlinear  solution  for  3-D  flows,  which 
is  expected  to  differ  only  slightly  from  the  fully  nonlinear  solution. 

One  virtue  of  using  only  the  first  iteration  solution  is  the  significant  computational 
savings  it  represents.  Both  the  fully  wetted  and  the  cavity  solutions  are  found  for 
a  fixed  discretization,  so  that  the  influence  coefficients  need  only  be  computed  once. 
The  resulting  3-D  solution  has  been  shown  to  capture  the  nonlinear  thickness  effect. 
The  dynamic  boundary  condition  has  been  verified  by  the  pressure  validation  test. 
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An  additional  time-saving  innovation,  which  also  solved  the  problem  of  poor  plan- 
form  convergence,  was  the  split-panel  method.  This  allowed  cavity  planforms  to  be 
continuous  and  general  without  adding  the  expense  of  computing  additional  influence 
coefficients. 

The  method  developed  in  this  thesis  for  treating  general  cavity  planforms  may  be 
described  as  a  uniform  formulation ,  meaning  that  partial  and  supercavity  planforms 
are  both  treated  within  the  same  mathematical  framework  using  the  same  discretiza¬ 
tion.  This  makes  the  treatment  of  mixed  cavity  planforms,  which  are  prevalent  on 
cavitating  propellers,  simple  and  natural. 

The  method  revealed  the  existence  of  multiple  solutions  for  high  aspect-ratio  rect¬ 
angular  hydrofoils.  This  multiplicity  is  reminiscent  of  the  analogous  2-D  case,  where 
the  flat  plate  solution  shows  that  there  are  three  cavity  lengths  for  each  cavitation 
number  for  ^  <  0.1.  For  this  case,  experiments  have  shown  that  2-D  cavities  longer 
than  l-  as  0.75  are  violently  unstable,  which  effectively  eliminates  the  occurrence  two 
of  the  three  predicted  cavity  lengths.  No  such  observations  have  been  reported  for 
3-D  hydrofoils.  Even  though  it  has  generally  been  believed  that  the  2-D  instability 
(and  corresponding  multiplicity)  disappear  due  to  “three  dimensional  effects”,  it  is 
highly  plausible  that,  for  high  aspect  ratio  3-D  foils,  neither  the  instability  nor  the 
multiplicity  of  solutions  are  mitigated  by  three-dimensionality.  The  multiplicity  has 
been  shown  to  disappear  as  the  aspect  ratio  is  lowered.  No  multiplicity  has  been 
found  for  the  propeller  solution,  although  a  systematic  investigation  has  yet  to  be 
tried. 

In  comparing  the  present  method  to  the  results  of  linear  theory  and  linear  theory 
with  leading  edge  corrections  for  two  dimensional  flows,  I  found  that  linear  theory 
substantially  overpredicted  the  cavity  length  and  volume  (this  was  not  a  new  result). 
The  linear  theory  with  leading  edge  corrections  was  much  closer,  but  at  high  angles 
of  attack  the  computed  cavity  volume  differed  substantially  from  the  exact  volume. 
Meanwhile,  even  for  high  angles  of  attack  the  first  iteration  of  the  nonlinear  solu¬ 
tion  compares  very  well  to  the  exact  solution.  As  a  result,  the  first  iteration  of  the 
nonlinear  solution  was  adopted  in  the  present  method. 
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In  making  a  similar  comparison  for  the  propeller  solution,  I  found  that  cavity 
planforms  predicted  by  linear  theory  and  linear  theory  with  leading  edge  corrections 
each  differed  from  the  nonlinear  result.  The  linear  theory  with  leading  edge  correc¬ 
tions,  while  underpredicting  the  cavity  extent,  seemed  to  capture  the  correct  character 
of  the  solution.  Meanwhile,  the  conventional  linear  theory  overpredicted  the  cavity 
planform  in  both  spanwise  and  chordwise  directions. 

The  pressure  validation  test  was  used  as  means  of  validating  the  results  of  the 
present  method  with  respect  to  the  fully  nonlinear  results.  The  pressure  validation 
test  showed  that  the  first  iteration  cavity  shapes  satisfy  the  dynamic  boundary  condi¬ 
tion  within  acceptable  accuracy,  from  which  it  is  concluded  that  the  shapes  are  close 
to  the  exact  nonlinear  shapes.  The  present  method  is  able  to  predict  cavity  shapes 
and  pressure  distributions  near  the  leading  edge  and  the  tip  with  greater  accuracy 
than  linear  theory,  especially  for  propellers  with  extreme  geometry.  We  conclude  that 
the  method  developed  in  this  thesis  is  preferable  to  that  of  PUF-3A.  However,  the 
present  method  labors  under  extensive  CPU  times  compared  to  PUF-3A.  For  example, 
a  fully  unsteady  solution  (through  120  time  steps,  or  two  revolutions)  of  a  one-bladed 
propeller  with  (N,M)=(80,16)  took  8  hours  on  a  DEC  9000.  PUF-3A,  which  uses  a 
discretization  of  40  chordwise  and  9  spanwise  discrete  vortex/source  lines,  took  2.5 
hours  for  the  same  geometry.  For  the  near  future,  the  present  method  is  more  likely 
to  be  used  as  a  validation  tool  for  PUF-3A.  However,' with  ever-improving  computer 
hardware  and  increasing  computer  speeds,  it  may  not  be  long  before  the  use  of  the 
method  as  a  design  tool  is  feasible.  One  need  only  look  at  the  CPU-time  history  of 
PUF-3A  as  a  recent  example  of  a  prohibitively  expensive  tool  becoming  cost-efficient. 


6.2  Recommendations  for  Future  Research 

The  following  improvements  may  be  made  to  the  current  method. 

•  The  method  is  in  need  of  an  efficient  algorithm  for  applying  the  nonlinear  pres¬ 
sure  Kutta  condition.  It  should  be  possible  to  implement  a  base-problem  solu¬ 
tion,  such  as  the  one  developed  by  Kinnas  and  Hsin  [38]. 
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•  The  first  dipole  panel  in  the  wake  should  be  given  a  linear  strength  distribution 
in  order  to  make  the  solution  more  nearly  independent  of  time  step  size  [22]. 
This  may  be  readily  implemented  in  the  present  code. 

•  The  ability  to  predict  face  cavitation,  where  the  supercavity  detaches  forward 
of  the  trailing  edge,  may  be  added  to  the  model.  It  may  also  be  possible  to 
treat  partial  cavities  which  occur  near  the  leading  edge  on  the  pressure  side. 
However,  this  would  be  logically  more  difficult  to  implement. 

Other  improvements  which  are  recommended  for  further  research  are  more  funda¬ 
mental  in  nature.  As  mentioned  in  Chapter  1,  the  fundamental  assumptions  made 
in  order  to  admit  a  potential  flow  solution  render  the  model  incomplete.  The  effects 
of  viscosity,  tip  and  hub  vortex  cavitation,  bubble  and  cloud  cavition,  and  the  cavity 
trailing  edge  flow  are  all  unaccounted  for  in  the  present  model.  However,  the  present 
model  is  amenable  to  inclusion  of  many  of  these  effects.  For  instance,  the  effects  of 
viscosity  may  be  included  via  an  interactive  viscous/inviscid  boundary  layer  solution, 
similar  to  the  one  developed  by  Hufford  for  fully  wetted  flows  [23].  The  boundary 
layer  solution  could  be  used  to  determine  the  thickness  of  the  cavity  wake  so  that  the 
openness  of  the  cavity  correctly  correlates  to  the  sectional  drag  coefficient.  A  model 
for  determining  the  correct  detachment  point,  by  correlating  the  point  of  laminar 
separation  or  turbulent  transition  to  the  point  of  cavity  detachment  (as  suggested 
by  several  previous  researchers  [2,  13]),  may  also  be  implemented.  A  model  of  the 
cavitating  tip  vortex  may  be  added  to  improve  the  solution  at  the  tip.  This  could  be 
accomplished  by  treating  the  tip  vortex  as  an  inner  problem,  the  solution  of  which 
should  be  matched  to  the  outer  solution  from  the  present  method.  The  inner  problem 
could  be  treated  by  a  boundary  element  method  with  a  grid  which  is  chosen  to  fit  a 
vortex  with  an  assumed  core  radius  (possibly  determined  semi-empirically). 

Before  implementing  any  of  these  models,  however,  a  systematic  series  of  experi¬ 
ments  should  be  performed  to  better  define  the  physical  nature  of  the  flow.  For  exam¬ 
ple,  detailed  flow  mapping  should  be  done  to  determine  the  nature  of  the  boundary 
layer  in  front  of  and  behind  the  cavity.  This  should  be  done  first  for  a  two-dimensional 
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hydrofoil.  It  would  also  be  useful  to  obtain  more  detailed  velocity  profiles  above  and 
approaching  the  cavity  surface,  especially  near  the  cavity  trailing  edge.  Experiments 
should  also  be  used  to  systematically  validate  the  computations  for  propeller  sheet 
cavitation.  Accurate  measurements  of  cavity  planforms  should  be  obtained,  first  for  a 
propeller  in  steady  flow,  for  comparison  to  computed  planforms.  Thrust,  torque,  effi¬ 
ciency  and  spanwise  circulation  distributions  may  also  be  measured  and  compared  to 
computed  values.  Later,  a  wake  screen  should  be  used  to  create  a  circumferentially 
nonuniform  inflow  to  the  propeller  and  a  similar  systematic  comparison  should  be 
made.  While  the  propeller  is  being  studied,  velocity  measurements  should  be  made 
through  a  tip  vortex  to  determine  its  size.  The  effect  of  parameters,  such  as  advance 
coefficient,  cavitation  number,  and  Reynolds  number  on  the  vortex  core  diameter 
should  be  investigated.  The  boundary  layer  near  the  tip  and  its  relationship  to  the 
vortex  core  size  should  also  be  studied. 

A  firm  foundation,  in  the  form  of  an  accurate  inviscid  solution,  is  required  to  build 
our  more  physical  models.  The  vortex/source  lattice  method,  as  implemented  in  the 
code  PUF-3A,  is  not  a  very  strong  foundation  due  to  the  neglect  of  blade  thickness 
and  the  crude  treatment  of  the  propeller  tip  geometry.  The  present  method,  however, 
is  an  accurate  inviscid  solution  and  will  serve  as  a  good  base  for  implementing  the 
additional  models. 
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Appendix  A 

Green’s  Formula  on  the  Wake  Surface 


The  goal  of  this  Appendix  is  to  prove  that  equation  (2.9)  of  Chapter  2  is  correct. 
This  equation  is  a  modified  version  of  Green’s  formula,  when  the  control  point  lies  on 
the  upper  half  of  the  wake  surface  (see  section  2.2).  The  canonical  problem  of  a  fully 
wetted  hydrofoil  with  a  wake  will  be  considered  (see  Figure  A-l).  The  wake  will  be 
a  surface  of  discontinuity  in  <f>  only,  since  the  treatment  of  the  source  discontinuity  in 
the  cavity  flow  problem  is  relatively  straightforward. 

Consider  a  hydrofoil  subject  to  an  incoming  free  stream  U-m,  as  shown  in  Figure 
A-l.  If  there  exists  a  velocity  potential  <j>  which  is  harmonic  everywhere  in  the  flow 
field,  then  it  is  a  classic  result  that  <f>  must  satisfy  Green’s  3r<*  identity,  which  can  be 
written  in  the  following  form 
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Figure  A-l:  Hydrofoil  and  wake,  showing  the  branch  cut  on  the  wake  surface. 
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Figure  A-2:  A  blow-up  of  the  field  point  p,  showing  the  small  patch  of  the  surface 
replaced  by  a  hemisphere  of  radius  R. 

where  p  and  q  are  the  field  point  and  the  variable  point  in  the  integration,  respectively. 
G(p;q)  =  1  /R(p;q)  is  Green’s  function,  where  R{p\q)  is  the  distance  between  the 
points  p  and  q.  The  surface  Sw  is  the  upper  half  of  the  wake,  as  shown  in  Figure 
A-l.  The  value  of  the  integral  I  depends  on  the  location  of  the  field  point: 

’ 

4ir<j>p  V  p  outside  Sr,  Sw 

I  =  '  0  V  p  inside  Sb  (A.2) 

2ir<f>p  V  p  on  Sb 


Equation  (A.l)  is  a  Fredholm  integral  equation  of  the  second  kind  when  the  field 
point  is  either  in  the  flow  field  or  on  the  hydrofoil  surface.  However,  equation  (A.2) 
does  not  define  I  for  the  case  when  the  field  point  lies  on  a  branch  cut,  which  is  the 
case  for  the  cavity  flow  problem  discussed  in  Chapter  2.  In  the  following,  equation 
(A.2)  will  be  generalized  to  include  the  case  where  p  lies  on  the  upper  side  of  a  surface 
across  which  4>  is  discontinuous. 

Consider  Green’s  3rd  identity  when  the  field  point  lies  on  the  wake  surface  Sw- 
In  this  case,  the  integral  over  Sw  must  exclude  the  field  point  in  order  for  Green’s 
formula  to  hold.  This  is  done  by  surrounding  the  field  point  by  a  hemisphere  Sr  of 
radius  R  which  indents  the  surface  into  the  external  flow,  as  shown  in  Figure  A-2. 
Green’s  formula  now  takes  the  form 


9-^r  -«(*»§]  dS+h  /  w-<>^r"=0'  <A-3> 

*  Sw+Sr  * 
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(A.4) 


Consider  the  integral  on  the  right-hand-side  of  equation  (A.4).  We  first  split  the 
integral  into  two  integrals 

5  /«  -  *;>T5T'S  ■  5  -  r.  /  (« 

5*  q  SR  q  SR  q 

In  the  limit  as  R  tends  to  zero  the  two  integrals  on  the  right-hand-side  of  (A. 5)  may 
be  evaluated  under  the  assumption  that  <f>  is  regular  on  Sr.  Equation  (A.5)  becomes 


—  f{<j>+  -  <f,~)?2kiSldS  =  --  lim 

4ir  J  'Vq  Vq)  dna  4tt  k-o 
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1  r 
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4ir  fl— o  dr  r 
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•2*  R2 
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•  2  *J*2  =  -±[#  +  <f>~] 


(A.6) 


J  r=fl 


Green’s  formula  for  the  case  where  the  field  point  lies  on  the  upper  wake  surface 
may  now  be  written: 


dS 


Sw  q 


Equation  (A. 2)  may  now  be  generalized  to 


/  = 


4  *<t>P 
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2tt^p 


V  p  outside  Sb,Sw 

V  p  inside  Sr 

V  p  on  Sr 


(A.7) 


2tt(0+  +  <f>~ )  V  p  on  Sw 
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Appendix  B 


Performance  of  the  Method  for  2-D 
Flows 


The  numerical  method  was  developed  first  in  two- dimensions  before  extending  it  to 
three  dimensional  geometries.  The  reasons  for  this  chronology  are  twofold:  first, 
it  is  useful  to  establish  convergence  and  consistency  of  the  method  in  2-D,  where 
algorithms  are  simple  and  computing  times  are  negligible;  and  second,  the  general 
characteristics  of  a  method  are  the  same  for  2-D  and  3-D  applications,  a  fact  which 
allows  the  basics  of  a  method  to  be  established  in  2-D  before  tackling  the  geometrically 
more  complex  3-D  problem.  For  example,  one  of  the  main  characteristics  of  the 
present  nonlinear  method  —  which  I  will  demonstrate  in  this  Appendix  —  is  that 
the  first  iteration  solution  is  very  close  to  the  converged  nonlinear  solution  for  a  wide 
range  of  operating  conditions.  This  may  be  said  to  be  a  characteristic  of  the  potential- 
based  panel  method,  since  the  velocity- based  panel  method  (whose  influence  functions 
are  one  degree  more  singular)  converges  very  slowly  with  iterations  [70].  Clearly,  the 
quick  convergence  with  iterations  is  a  characteristic  which  may  be  expected  to  hold 
for  3-D  geometries. 

In  addition  to  demonstrating  the  above-mentioned  rapid  convergence,  the  objec¬ 
tive  of  this  Appendix  is  to  display  numerical  convergence  and  consistency  of  the  2-D 
solution.  Details  of  the  solution,  which  are  entirely  analogous  to  the  details  of  the 
unsteady  propeller  solution  described  in  the  main  body  of  this  thesis,  are  reported  in 
several  articles  which  are  included  in  the  references. 
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nr 

<7 

mm 

0.0036 

H 

0.4493 

0.0034  | 

El 

0.4502 

0.0033 

100 

0.0033 

120 

0.0033 

Table  B.l:  Convergence  with  number  of  panels  for  a  partially  cavitating  NACA16006 
symmetric  foil  at  a  =  3°,  £  =  0.5  and  l-f  =  0.048.  The  termination  model  parameters 
are:  j  =  0.1,  A  =  0.3  and  v  =  1.0  (for  the  definition  of  these  parameters,  see 
Appendix  C). 


B.l  Fixed-length  Solution 

The  investigation  of  the  numerical  method  began  with  the  two-dimensional  fixed- 
length  problem,  as  reported  in  [35].  In  what  follows,  I  will  show  the  results  of  conver¬ 
gence  (with  number  of  panels  and  with  iterations)  and  consistency  tests  for  partially- 
and  super-cavitating  2-D  hydrofoils  with  fixed-length  cavities. 


B.1.1  Partial  Cavitation 

Table  B.l  shows  convergence  with  number  of  panels  for  the  first  iteration  solution  (all 
panels  on  the  foil  surface)  of  a  partially  cavitating  (l  —  0.4)  NACA16006  symmetric 
hydrofoil  at  a  =  3°.  There  are  two  things  to  learn  from  this  convergence  test.  The 
first  and  most  important  is  simply  that  the  solution  does  converge.  The  second  is 
the  number  of  panels  required  to  obtain  what  I  consider  to  be  a  converged  solution 
(three  significant  figures  in  <r).  In  this  case  that  number  is  near  80,  which  is  more 
than  that  required  for  the  fully-wetted  solution  (which  requires  roughly  60  panels  for 
the  circulation  to  converge  to  3  significant  figures  [47]).  One  likely  explanation  for 
needing  more  panels  for  the  cavity  solution  is  the  sensitivity  to  the  leading  edge  flow 
[33,  70].  I  have  found  that  clustering  panels  near  the  leading  edge  (at  the  expense  of 
panel  density  aft  of  the  leading  edge)  sped  the  convergence  of  the  fixed  length  solution 
slightly. 

Figure  B-l  shows  the  solution  history  as  it  proceeds  from  the  first  iteration  solution 
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Figure  B-l:  Convergence  history  of  a  fully  nonlinear  partial  cavitation  fixed-length 
solution.  NACA16006  foil  at  a  =  3°,  j  =  0.4  and  k  =  0.0.  The  termination  model 
parameters  are  A  =  0.1,  A  —  0.5,  and  v  —  1.0;  the  number  of  panels  is  N  =  100.  The 
cavitation  number  and  cavity  volume  shown  in  the  plot  on  the  right  are  normalized 
on  their  converged  values. 

to  the  fully  nonlinear  solution.  The  progression  of  cavity  shapes  is  shown  on  the  left 
of  the  figure  and  the  cavitation  numbers  and  the  cavity  volumes  (normalized  on  the 
converged  values)  are  shown  on  the  right.  The  exceptionally  quick  convergence  which 
is  depicted  in  this  figure  is  characteristic  of  the  solution  at  a  wide  range  of  angles  of 
attack  and  foil  thicknesses.  This  will  be  shown  for  the  solution  of  the  fixed-<r  problem, 
which  is  described  in  the  next  section. 

Figure  B-2  shows  pressure  validation  of  the  first  and  last  iterations  in  the  nonlinear 
solution  for  a  NACA16006  foil  at  a  =  3°  and  ~  =  0.4.  The  pressure  validation,  which 
is  also  described  in  section  5.2.2,  consists  of  computing  the  fully  wetted  flow  about 
a  modified  foil  made  up  of  the  original  foil  and  the  computed  cavity.  The  resulting 
pressure  distribution  should  be  constant  and  equal  to  —  <r  on  the  cavity.  The  deviation 
from  this  is  a  measure  of  the  error  introduced  by  iterating  only  once.  Notice  that  the 
pressures  match  exactly  for  the  fully  nonlinear  solution  and  that,  even  for  the  first 
iteration  solution,  the  pressures  sure  very  close. 
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Figure  B-2:  Pressure  validation  test  for  the  first  (top)  and  last  (bottom)  iterations 
of  the  nonlinear  solution  presented  in  the  previous  figure.  Shown  to  the  left  of  each 
pressure  distribution  is  the  corresponding  panel  arrangment  for  the  fully  wetted  so- 


N 

a 

Volume/ c2 

80 

0.130 

0.190 

120 

0.124 

0.183 

160 

0.122 

0.180 

200 

0.123 

0.181  | 

240 

0.123 

0.181 

Table  B.2:  Convergence  with  number  of  panels  of  the  fixed-length  solution  for  a 
supercavitating  NACA16004  symmetric  foil  at  a  =  5°,  1  =  2.0  and  k  =  0.024.  The 
termination  model  parameters  are:  =  0.07,  ^  =  0.07  (for  the  definition  of  these 

parameters,  see  Appendix  C).  The  number  of  panels  is  N  —  100. 

B.1.2  Supercavitation 

The  nonlinear  fixed-length  solution  for  supercavitating  hydrofoils  is  found  by  using 
the  cavity  shape  predicted  by  linear  theory  as  the  initial  guess.  This  is  a  much  better 
initial  guess  than  the  one  used  for  the  partial  cavitation  fixed -l  solution,  since  linear 
theory  is  known  to  be  more  accurate  for  long  cavities  than  for  short.  Note  that 
now  the  phrase  “first  iteration”  means  that  the  panels  are  arranged  on  the  surface 
determined  from  linear  theory  (so  the  linear  theory  is  the  zeroth  iteration).  Even 
though  there  is  no  reasonable  3-D  extension  to  an  algorithm  which  uses  linear  theory 
as  an  initial  guess,  it  is  nonetheless  useful  to  investigate  the  character  of  this  solution. 
For  example,  the  fixed-length  solution  will  be  used  to  investigate  termination  models 
in  Appendix  C.  Meanwhile,  since  the  code  is  available,  it  does  not  hurt  to  verify  that 
the  solution  is  convergent  and  consistent. 

Table  B.2  shows  convergence  with  equivalent  number  of  panels.  N  is  the  number 
of  panels  distributed  around  the  entire  foil/cavity  surface,  where  the  cavity  surface 
is  the  one  predicted  by  linear  theory.  The  cavity  length  is  £  =  2.0,  so  the  equivalent 
number  of  panels  is  twice  that  in  Table  B.l. 

Figure  B-3  shows  convergence  with  iterations  of  the  cavity  shape,  cavitation  num¬ 
ber  and  cavity  volume  for  a  supercavitating  NACA16006  foil  at  a  =  3°,  ~  =  1.5 
and  k  =  0.05.  In  this  plot,  the  cavity  shapes  are  shown  from  linear  theory,  the 
nonlinear  1**  iteration,  and  the  converged  nonlinear  solution.  The  shapes  from  all 
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Figure  B-3:  Convergence  of  cavity  shapes,  cavitation  number,  and  cavity  volume 
with  number  of  iterations  for  a  supercavitating  NACA16006  foil  at  a  =  3°,  ~  =  1.5, 
k  =  0.05,  =  0.05,  ^  =  0.05.  The  number  of  panels  is  N  =  100. 

three  solutions  are  nearly  indistinguishable.  Figure  B-4  shows  pressure  validation  for 
the  foil  and  conditions  shown  in  Figure  B-3.  Here,  only  the  first  iteration  pressure 
validation  is  shown,  since  the  other  solutions  are  nearly  identical.  Note  that  the 
pressure  distribution  reveals  pressures  on  the  lower  surface  which  are  lower  than  the 
cavity  pressure.  This  means  that  part  of  the  lower  surface  near  the  foil  trailing  edge 
should  be  cavitating  and  underscores  the  need  to  include  so-called  “face  cavitation” 
in  the  model.  This  was  done  using  linear  theory  by  Kinnas  and  Fine  [38]  (where  the 
suction-side  cavity  detachment  point  was  found  from  the  Villat-Brillouin  condition), 
and  should  be  added  to  the  present  method.  Figure  B-5  shows  convergence  with 
iterations  for  the  same  foil  and  conditions  as  in  Figure  B-3,  except  that  in  this  case, 
the  cavity  detaches  at  the  leading  edge  (^  =  0.0).  Note  that  the  three  solutions  are 
now  different  and  the  linear  theory  no  longer  comes  close  to  predicting  the  correct 
cavity  shape.  This  may  be  due  to  the  fact  that  the  exact  cavity  surface  intersects  the 
foil  near  the  leading  edge,  as  shown  in  Figure  B-5. 
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Figure  B-4:  Pressure  validation  of  a  supercavitating  NACA16006  foil.  Same  foil  and 
conditions  as  in  Figure  B-3. 

B.2  Fixed-cr  Solution 

In  practical  problems,  including  the  problem  of  predicting  unsteady  propeller  sheet 
cavitation,  the  cavitation  number  is  known  and  the  cavity  length  is  to  be  determined. 
In  light  of  the  ultimate  goals  of  the  research,  it  was  clearly  beneficial  to  solve  this  so- 
called  “direct  problem”  in  two  dimensions  before  attempting  it  in  three.  The  details 
of  the  algorithm  for  the  first  iteration  solution  have  been  described  in  the  main  body 
of  this  thesis  and  in  [37]  and  [13].  The  details  of  the  fully  nonlinear  solution  will 
be  given  in  section  B.2.3.  In  this  section,  I  will  show  results  of  convergence  and 
consistency  tests  for  the  2-D  fixed-<7  solutions. 

B.2.1  Partial  Cavitation 

In  all  of  the  convergence  and  consistency  tests  shown  in  this  section,  the  solution 
is  found  for  constant  panelling  (using  the  split  panel  method).  Table  B.3  shows 
convergence  of  cavity  length  and  volume  with  number  of  panels  (first  iteration  solution 
only)  for  a  NACA16006  symmetric  foil  at  a  =  3°.  There  are  two  reasons  why  the 
fixed- <t  solution  does  not  converge  uniformly  with  number  of  panels:  first,  there  is  a 
small  amount  of  error  due  to  the  split  panel  method,  and  that  error  does  not  varnish 
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Figure  B-5:  Convergence  of  cavity  shapes  with  number  of  iterations  for  the  same 
foil  and  conditions  as  the  previous  figure,  except  that  k  =  0.0.  In  the  blow  up  of 
the  leading  edge,  it  is  clear  that  the  nonlinear  solutions  cut  the  foil,  while  the  linear 
solution  does  not. 
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N 

l 

C 

Volume/c2 

■ 

1 

mM 

0.3852 

0.00366 

4.0  x  10-6 

80 

0.3944 

0.00403 

7.9  x  10-5 

100 

0.3915 

0.00392 

8.3  x  lO"5 

120 

0.00409 

5.8  x  lO"5 

140 

0.3954 

0.00400 

-1.9  x  lO'5 

Table  B.3:  Convergence  with  number  of  panels  of  the  partial  cavitation  fixed- 
(Tsolution  for  a  NACA16006  symmetric  foil  at  a  =  4°,  a  =  0.7  and  ^  =  0.024. 
The  termination  model  parameters  are:  j  =  0.1,  A  =  0.3,  and  u  =  1.0.  The  panels 
are  arranged  with  cosine  spacing,  as  shown  above  for  N=60. 

locally  with  increased  number  of  panels.  Thus  a  (very)  slightly  different  problem  is 
being  solved  every  time,  making  the  convergence  appear  slow  or  nonuniform.  Second, 
the  solution  is  found  to  within  a  tolerance  on  6 ,  and  the  final  tolerance  may  vary  from 
one  solution  to  the  next.  Invariably,  the  coarser  the  grid,  the  larger  the  tolerance  must 
be  in  order  to  arrive  at  a  solution.  Despite  these  handicaps,  Table  B.3  shows  that  the 
method  gives  the  same  cavity  length  and  volume  (±2.5%)  for  all  N  greater  than  or 
equal  to  80.  The  error  caused  by  the  split  panel  method  has  been  discussed  in  section 
4.1.1. 

Another  important  test  for  the  fixed-<r  solution  is  to  vary  the  initial  guess  and 
check  that  the  final  cavity  lengths  are  close  to  one  another.  This  test  is  a  gauge  of  the 
error  introduced  by  split  panel  method.  Table  B.4  shows  such  a  test  for  the  same  foil 
tested  in  Table  B.3.  Shown  also  is  the  correct  cavity  length  (which  is  independent  of 
initial  guess),  found  by  treating  the  split  panel  as  two  panels. 

Figure  B-6  shows  the  first  and  last  iterations  in  convergence  to  the  fully  nonlinear 
solution  (see  section  B.2.3  for  details).  Also  shown  is  the  corresponding  pressure 
validation  from  each  solution.  This  figure  shows  that  the  first  iteration  cavity  solution 
is  very  close  to  the  converged  nonlinear  result. 

Figure  B-7  shows  a  comparison  between  the  results  of  linear  theory,  linear  the¬ 
ory  with  leading  edge  corrections,  the  present  method  (first  iteration)  and  the  fully 
nonlinear  solution.  Shown  are  cavity  shapes  from  each  solution,  superimposed  on 
one-another.  The  linear  theory  is  seen  to  overpredict  the  cavity  extent  and  volume 
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Figure  B-6:  Partially  cavitating  NACA16006  foil  at  a  =  3°,  A  =  0.1,  v  —  1.0,  A  —  0.5, 
or  =  .731.  The  top  plot  shows  the  first  and  last  iterations  in  the  nonlinear  fixed-<7 
solution.  The  last  iteration  corresponds  to  a  cr-tolerance  of  1%.  Also  shown  is  the 
converged  nonlinear  ftxed-f  solution.  The  bottom  plot  shows  pressure  validation  for 
the  first  iteration  and  the  fully  nonlinear  solutions. 
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Table  B.4:  Varying  the  initial  guess  for  the  same  foil  and  conditions  as  in  the  Table 
A.3.  The  number  of  panels  is  N=80,  and  the  panel  arrangement  is  shown  above. 


for  each  operating  condition.  The  linear  theory  with  leading  edge  corrections  predicts 
the  cavity  extent  well,  but  overpredicts  the  volume  for  high  angles  of  attack  and  small 
thicknesses.  On  the  other  hand,  the  first  iteration  nonlinear  shapes  are  very  close  to 
the  converged  nonlinear  shapes  for  each  angle  of  attack  and  each  foil  thickness. 
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Figure  B-T:  Comparison  of  cavity  shapes  predicted  from  different  methods  for  a 
NACA16  foil  at  a  —  4°  (top),  6°  (middle),  and  8°  (bottom),  and  ~  =  .06  (left)  and 
.09  (right).  The  ratio  f  is  kept  fixed  at  .055.  The  cavity  detaches  at  the  leading  edge 
in  all  cases. 


B.2.2  Supercavitation 

The  fixed-<r  supercavitating  hydrofoil  problem  has  been  solved  only  for  the  first  it¬ 
eration,  where  Green’s  formula  is  satisfied  on  the  foil  surface  and  on  the  portion  of 
the  wake  which  is  overlapped  by  the  supercavity.  Finding  the  fully  nonlinear  solution 
in  this  case  was  deemed  unnecessary  because  nothing  would  be  learned  about  the 
method  that  could  not  be  inferred  from  the  partial  cavitation  fixed-<r  or  the  super¬ 
cavitating  fixed-length  solutions. 

Table  B.5  displays  convergence  with  number  of  panels  for  a  supercavitating  sym¬ 
metric  NACA16006  foil  at  a  =  6°.  There  are  N  panels  distributed  in  the  usual 
cosine-spacing  on  the  foil  and  Nc  constant-spaced  panels  in  the  wake.  The  wake  is 
re-panelled  after  each  guess  of  the  cavity  planform,  so  that  the  cavity  always  ends  at 
a  panel  boundary  ( e.g .,  there  is  no  split  panel). 

Figure  B-8  shows  a  comparison  of  cavity  solutions  from  linear  theory,  the  present 
method  first  iteration  and  the  present  method  converged  nonlinear.  Here  the  con¬ 
verged  nonlinear  shape  is  found  from  the  fixed-length  solution,  and  the  other  solutions 
are  found  for  the  a  which  resulted  from  the  fixed-/  nonlinear  solution.  Cavity  shapes 
are  shown  for  a  NACA16006  foil  at  a  =  4°,  6°,  and  8°. 

B.2.3  Partial  Cavitation  Fully  Nonlinear  Solution 

It  is  not  immediately  clear  how  to  proceed  from  the  first  iteration  cavity  length  to  a 
converged  nonlinear  cavity  length  for  the  fixed- cr  solution.  If  we  “build  on”  the  first 
iteration  cavity  by  moving  the  panels  which  represent  the  cavity  from  the  foil  surface 
to  the  newly  computed  cavity  surface  and  apply  Green’s  formula  there  for  the  given 
cavitation  number,  then  the  result  will  be  a  cavity  which  does  not  close.  In  fact, 
since  we  know  from  experience  with  the  fixed-length  solution  that  the  nonlinear  <r  is 
always  less  than  the  first  iteration  <t,  we  can  predict  that  the  first  iteration  cavity  is 
too  long.  Therefore,  building  on  the  first  iteration  cavity  usually  results  in  a  cavity 
for  which  S  <  0.  By  trimming  off  the  end  of  the  cavity  which  lies  inside  the  foil,  the 
cavity  length  may  be  shortened  and  another  iteration  may  be  built  on  the  resulting 
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N 

Nc 

l 

Volume/c2 

& final 

KTj] 

El 

1.6348 

■ 

4.3  x  10-6 

■Si 

El 

1.6285 

3.8  x  10-* 

60 

H 

1.6262 

0.1391 

3.3  x  10“8 

80 

10 

0.1415 

4.6  x  10"8 

80 

20 

t‘  ___ 

0.1408 

3.5  x  10-8 

80 

40 

3.5  x  10-8 

100 

10 

1.6415 

0.1415 

4.7  x  10-8 

100 

20 

1.6348 

0.1407 

3.8  x  10-8 

100 

40 

1.6323 

0.1403 

3.4  x  10-8 

120 

10 

0.1412 

4.7  x  10-8 

120 

20 

1.6343 

0.1404 

3.8  x  10-8 

120 

40 

1.6317 

0.1401 

3.5  x  10'8 

120 

60 

1.6310 

0.1400 

3.4  x  10-8 

Table  B.5:  Convergence  with  number  of  panels  for  a  NACA16006  symmetric  foil  at 
a  =  6°,  a  —  0.2  and  k  =  0.0. 
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Figure  B-8:  Comparison  of  supercavity  shapes  predicted  from  linear  theory,  the 
present  method,  and  nonlinear  theory  for  a  NACA16004  foil  at  a  =  4°  (top), a  =  6° 
(middle),  and  a  =  8°  (bottom).  The  ratio  f  is  kept  fixed  at  0.440.  In  each  case,  the 
cavity  detaches  at  the  leading  edge. 
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Figure  B-9:  An  example  of  the  original  method  for  finding  the  nonlinear  fixed-er 
solution.  The  convergence  is  seen  here  to  be  very  slow. 

cavity  shape.  Repeating  this  procedure  will  eventually  lead  to  the  nonlinear  cavity 
shape.  Unfortunately,  I  found  that  the  convergence  of  this  scheme  was  very  slow,  as 
is  evident  in  Figure  B-9. 

An  accelerated  scheme  has  been  developed  which  takes  advantage  of  the  fact 
that  the  difference  between  the  cavitation  numbers  from  the  first  and  the  converged 
iterations  with  fixed  cavity  length,  |,  behaves  smoothly  with  With  or  being  the 
given  cavitation  number,  the  proposed  algorithm  is  as  follows: 

1.  Solve  the  cavity  problem  (first  iteration)  with  fixed  cavitation  number 
(Tq  —  (T  to  find  the  corresponding  cavity  length  l0. 

2.  Solve  the  cavity  problem  with  fixed  cavity  length  l„  for  K  iterations  (t.e. 
regridding)  and  find  the  converged  value  of  the  cavitation  number,  <r'0. 

3.  Define  a  new  cavitation  number,  <T\  =  <7o  4-  (<To  —  ^0)' 

4.  Solve  the  cavity  problem  (first  iteration)  with  fixed  cavitation  number  <rn 
and  find  the  correct  cavity  length  ln  ( n  >  1). 

5.  Solve  the  cavity  problem  with  fixed  cavity  length  ln  for  K  iterations,  and 
find  the  converged  value  of  the  cavitation  number,  <r'n  (n  >  1). 

6.  Apply  a  Newton-Raphson  (secant)  iterative  solution  for  <r'n  —  a  —  0,  and 
find  a  new  cavitation  number,  <rn+i . 
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7.  Repeat  steps  4,  5,  and  6  until  <r'n  =  <r  to  within  a  preset  tolerance. 


For  the  results  shown  in  Figure  B-6  the  last  iteration  was  found  for  a  tolerance  of 
2*^ -  <  .01.  Convergence  to  this  tolerance  took  only  one  Newton-Raphson  iteration. 
Since  the  nonlinear  cavitation  number  at  each  cavity  length  was  found  by  iterating 
twice  ( K  =  2),  the  total  number  of  regriddings  for  this  case  was  equal  to  4.  Included 
in  Figure  B-6  is  the  converged  nonlinear  cavity  (known,  since  a  —  0.731  was  chosen 
to  correspond  to  the  converged  cavity  solution  with  fixed  length  /  =  0.4).  Also  shown 
in  Figure  B-6  is  the  pressure  distribution  on  the  foil  and  cavity  from  the  converged 
nonlinear  solution  (plotted  with  a  solid  curve).  Superimposed  on  this  curve  are  two 
pressure  distributions  which  resulted  from  the  fully  wetted  analysis  of  the  foil  and 
cavity  from  the  first  iteration  (plotted  with  asterisks)  and  the  last  iteration  (plotted 
with  open  circles).  Note  that  the  pressure  distribution  on  the  cavity  from  the  last 
iteration  satisfies  the  imposed  dynamic  boundary  condition,  and  even  the  pressure 
distribution  on  the  cavity  from  the  first  iteration  comes  very  close  to  satisfying  the 
condition. 

B.3  Effect  of  the  detachment  point 

It  was  mentioned  in  Chapter  1  that  the  point  at  which  the  cavity  detaches  depends 
on  the  nature  of  the  viscous  boundary  layer.  Since  we  don’t  include  viscosity  in  the 
present  method,  the  point  of  detachment  is  considered  as  an  independent  variable  in 
the  solution.  In  this  section,  the  effect  of  the  detachment  point  will  be  investigated 
for  two  dimensional  flows. 

Figure  B-10  shows  the  first  iteration  solution  for  a  partially  cavitating  NACA16006 
foil  at  a  =  3°  with  several  points  of  detachment  for  fixed  cavitation  number  a  =  0.7. 
Here  it  is  seen  that  the  detachment  point  has  a  dramatic  effect  on  the  cavity  extent, 
with  the  cavity  length  diminishing  as  the  detachment  point  moves  aft. 
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Figure  B-10:  The  effect  of  detachment  point  on  the  partial  cavitation  solution.  Shown 
are  four  solutions  corresponding  to  detachment  points  ^  =  0,  .004,  .016  and  .035.  The 
foil  is  NACA16006  at  a  =  3°  and  a  =  0.7.  The  termination  model  parameter  A  is 
zero. 

B.4  Conclusions  from  2-D  Survey 

The  following  conclusions  may  be  drawn  from  my  survey  of  the  present  method  for 
2-D  flows: 

1.  The  nonlinear  cavity  shapes  from  the  first  and  converged  iterations  are 
very  close  to  one  another  for  a  wide  range  of  conditions.  This  suggests 
that  only  a  small  amount  of  error  is  introduced  by  satisfying  the  boundary 
conditions  on  the  approximate  flow  boundary.  The  cavity  shapes  predicted 
by  linear  theory  with  leading  edge  corrections  are  also  close  to  the  nonlinear 
shapes  for  small  angles  of  attack  and  large  thicknesses,  but  they  diverge 
from  the  nonlinear  shapes  as  the  angle  of  attack  increases  and  the  thickness 
decreases. 

2.  Typically  the  number  of  panels  required  for  3-digit  convergence  of  <r  is 
N  =  80. 

3.  An  additional  requirement  for  convergence  with  number  of  panels  is  the 
incorporation  of  an  accurate  definition  of  the  potential  at  the  leading  edge 
of  the  cavity,  fo. 
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Appendix  C 


Termination  Models 


The  need  for  a  termination  model  arises  from  the  inconsistency  inherent  in  forcing  a 
constant-velocity  streamline  to  end  at  a  stagnation  point.  This  need  has  long  been 
recognized.  As  mentioned  in  Chapter  1,  excellent  reviews  of  existing  models  may 
be  found  in  [75]  and  [63].  The  termination  models  investigated  for  this  thesis  are 
described  in  this  appendix. 

Many  termination  models  are  designed  first  and  foremost  for  mathematical  expe¬ 
diency.  Accurate  modelling  of  the  physics  tends  to  be  a  secondary  priority,  mainly 
because  it  is  impossible  to  simulate  the  physical  nature  of  multiphase  turbulent  flow 
using  a  potential  flow  model.  Velocity  measurements  made  behind  a  partial  cavity 
on  a  two-dimensional  foil  in  steady  flow  showed  the  existence  of  a  rather  thick  wake 
trailing  the  cavity  [11].  Measurements  on  the  cavity  also  showed  a  slight  velocity 
attenuation  near  the  cavity  trailing  edge.  These  experiments  suggest  that  an  open 
wake  model  may  be  the  most  physically  realistic  termination  model. 

In  this  thesis,  several  models  were  tried.  These  fall  into  two  categories: 

1.  The  end-plate  model.  This  model  is  similar  to  the  Riabouchinsky  end- plate 
model,  as  applied  by  Uhlman  [69],  except  that  the  end- plate  is  allowed  to 
have  curvature.  The  source  strength  or,  alternatively,  the  cavity  height, 
near  the  end  of  the  cavity  is  specified  by  an  algebraic  law. 

2.  The  pressure  model.  In  this  model,  first  suggested  by  Lemonnier  and  Rowe 
[20],  the  pressure  on  the  cavity  is  allowed  to  vary  in  a  region  near  the  end 
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Figure  C-l:  A  depiction  of  the  cavity  height  law  for  partial  cavitation. 

of  the  cavity  from  its  constant  value  upstream  to  some  fraction  of  that 
value  at  the  cavity  trailing  edge.  This  model  also  allows  for  the  cavity  to 
be  open  and  seems  to  be  a  more  physically  realistic  simulation. 

At  the  time  these  were  developed  for  2-D  flows,  the  motivation  was  to  find  a  model 
which  allowed  for  smooth  and  quick  convergence  to  the  fully  nonlinear  solution.  Some 
details  and  results  of  the  investigation  will  be  given  in  the  following  sections. 


C.X  Partial  Cavitation 


C.1.1  Cavity  Height  Law 

In  this  model,  the  cavity  is  split  into  two  zones:  one  zone  near  the  trailing  edge 
of  the  cavity  (called  the  transition  zone)  with  horizontal  projected  length  A,  with  the 
other  zone  being  the  upstream  remainder  (see  Figure  C-l).  Within  the  transition 
zone,  the  source  strength  is  specified  and  given  by  the  algebraic  expression 


9* ,  .  _  H .  , 


Si  —  s 


si  —  srJ 


st  <  s  <  si 


(C.l) 


Here,  s  is  the  arclength  of  the  foil  beneath  the  cavity,  measured  from  the  cavity  leading 
edge  along  the  foil  surface,  si  is  the  value  of  *  at  the  trailing  edge  of  the  cavity,  and 
st  is  the  value  of  s  at  the  beginning  of  the  transition  zone.  The  streamwise  extent  of 
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the  transition  zone  is  defined  by  the  parameter  A  =  (si,  —  sj )/&l-  The  parameters  A 
and  v  (v  >  0)  are  arbitrary  constants. 

An  example  of  a  computation  made  with  this  termination  model  is  shown  in 
Figure  C-2.  Shown  is  the  first  iteration  cavity  for  a  NACA16006  foil  at  a  =  4°  with 
a  cavity  length  of  l-  =  0.4.  The  length  of  the  transition  zone  is  £  =  0.15.  Although 
the  cavity  looks  smooth  and  reasonable,  the  method  suffered  from  poor  convergence 
to  the  fully  nonlinear  cavity  (which  is  why  I  only  show  the  first  iteration).  One 
possible  explanation  for  this  is  the  slope  discontinuity  of  the  cavity  at  s  —  st-  A 
slope  continuity  condition  could  have  been  incorporated.  However,  a  natural  way  to 
ensure  slope  continuity  is  to  use  a  model  which  governs  the  cavity  pressure  rather 
than  the  cavity  height.  The  height  model  was  then  abandoned  in  favor  of  a  pressure 
model,  which  will  be  described  next. 


C.1.2  Pressure  Law 


The  pressure  law  termination  model  was  developed  in  order  to  ensure  slope  conti¬ 
nuity  of  the  cavity  surface,  as  mentioned  above.  It  was  also  our  intention  to  simulate 
the  real  attenuation  of  pressure  which  has  been  measured  in  steady  flow  2-D  experi¬ 
ments  at  MIT  [24].  The  cavity  is  again  broken  into  two  zones,  with  a  transition  zone 
of  length  A  (see  Figure  C-3).  The  model  requires  the  total  velocity  on  the  cavity  to 
satisfy 

*  =  «.  [1  -  /«]  (C.2) 


where 


/(*)  = 


0  if  s  <  aj 


(C.3) 


The  parameters  are  defined  in  section  C.1.2,  with  the  exception  of  A,  which  is  an 
arbitrary  constant  between  zero  and  one.  A,  A  and  u  are  inputs. 

The  convergence  of  the  scheme  to  the  fully  nonlinear  solution,  with  the  pressure 
law  termination  model  implemented,  is  excellent.  This  has  already  been  shown  in 
Chapter  1  and  in  Appendix  B  (see,  for  example,  Figure  B-l). 
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Figure  C-2:  An  example  of  a  partial  cavity  solution  using  the  height  law  termination 
model.  Shown  are  the  cavity  shape  and  pressure  distribution  for  a  N  AC  A 16006  foil  at 
a  =  4°  with  a  cavity  length  of  ~  =  0.4.  The  length  of  the  transition  zone  is  £  =  0.15. 
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Figure  C-3:  A  depiction  of  the  pressure  law  for  partial  cavitation. 


Figure  C-4  shows  the  effect  of  varying  the  length  of  the  transition  zone,  A,  on  the 
converged  nonlinear  cavity  shape.  Varying  A  from  .06  to  .24,  with  all  other  parameters 
frozen,  decreased  the  volume  by  12%.  However,  for  the  same  perturbation  of  A,  the 
cavitation  number  increased  by  only  3%.  The  parameter  v  was  also  varied  in  a 
numerical  experiment,  which  showed  the  dependence  to  be  relatively  minor,  as  is  also 
evident  in  Figure  C-4. 

Figure  C-5  shows  the  dependence  of  the  solution  on  the  parameter  A.  Here,  it  is 
seen  that  the  solution  depends  strongly  on  A.  If  A  is  too  large  for  the  given  A,  then 
the  cavity  surface  intersects  the  foil  surface.  For  A  ~  0.15,  A  =  .5  usually  produces 
good  results. 


C.2  Supercavitation 

C.2.1  Curved  End-Plate  Model 

Several  end  plate  models  were  examined  in  the  course  of  the  investigation.  First,  a 
modified  Riabouchinsky  end  plate  model  was  tried,  similar  to  the  one  used  by  Uhlman 
[71].  In  this  model,  the  end  of  the  cavity  is  replaced  by  a  vertical  flat  plate  on  which 
the  kinematic  boundary  condition  is  satisfied.  The  height  of  the  plate  comes  out  of 
the  solution.  While  the  x  location  of  the  plate  was  fixed,  it  was  unconstrained  in 
the  vertical  plane.  Since  there  was  no  closure  condition  in  this  model,  an  additional 
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Figure  C-4:  Cavity  shapes,  cavitation  numbers  and  cavity  volumes  for  different  values 
of  the  parameters  A  and  v. 
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Figure  C-5:  Cavity  shapes  for  different  values  of  the  parameter  A.  NACA16006  foil 
with  q  =  4°  and  £  =  0.5. 
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Figure  C-6:  A  sketch  of  the  end-parabola  model  for  supercavitation. 

equation  was  needed  to  round  out  the  system  of  equations.  For  this  reason,  a  matching 
slope  condition  was  applied  at  the  trailing  edge  of  the  foil  (matching  the  pressure  side 
foil  slope  to  the  slope  of  the  cavity).  Unfortunately,  this  turned  out  to  be  an  ill 
posed  problem.  The  slope  matching  condition  was  unnecessary  because  linear  theory, 
which  was  used  as  the  first  iteration,  satisfies  a  Kutta  condition  at  the  foil  trailing 
edge  which  resulted  in  continuous  slopes.  Thus  the  problem  was  over-constrained 
and  the  resulting  cavity  shapes  were  unrealistic.  This  prompted  us  to  change  to  a 
model  which  used  a  different  kind  of  end-plate  model,  one  which  utilized  a  form  of 
the  cavity  closure  condition.  In  this  model,  the  end  plate  model  was  allowed  to  have 
curvature,  and  we  refer  to  it  as  the  curved  end  plate  model. 

The  curved  end  plate  model  (shown  in  Figure  C-6)  is  similar  to  the  cavity  height 
model  for  partial  cavitation.  This  model  uses  an  algebraically  defined  curved  end 
plate  on  which  the  kinematic  boundary  condition  is  applied.  The  end  plate  has 
horizontal  extent  defined  by  the  parameter  A,  as  shown  in  Figure  C-6.  In  this  case  a 
cavity  closure  condition  is  applied  which  constrains  the  end  plate  to  move  vertically. 
Specifically,  the  closure  condition  requires  the  upper  and  lower  halves  of  the  end  plate 
to  move  in  unison.  In  an  early  version,  the  curved  end  plate  was  simply  a  parabola 
which  maintained  its  shape  throughout  all  iterations;  it  simply  moved  up  or  down 
according  to  the  solution.  This  model  was  very  successful  in  obtaining  the  converged 
nonlinear  solution  in  relatively  few  iterations.  This  model,  however,  suffered  due  to 
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the  lack  of  a  requirement  that  the  slope  of  the  cavity  be  continuous  at  the  juncture 
of  the  end  parabola  and  the  upstream  part  of  the  cavity.  The  resulting  discontinuity 
in  slope  was  not  strong  enough  to  corrupt  the  iterative  convergence;  however,  it  did 
result  in  strange  looking  cavities  at  high  angles  of  attack.  Results  for  a  supercavitating 
NACA16006  foil  are  shown  in  Figure  C-7.  In  this  figure,  the  solution  is  shown  from 
linear  theory  and  nonlinear  theory  (using  the  end  parabola  model)  at  two  angles  of 
attack  (10°  and  18°).  At  the  higher  angle  of  attack,  the  discontinuity  is  clear. 

The  next  modification  was  to  include  matching  slope  conditions  to  ensure  slope 
continuity  at  the  cavity-parabola  juncture.  Rather  than  apply  the  conditions  as  two 
additional  equations,  which  would  have  rendered  the  problem  over-constrained,  we 
chose  to  satisfy  the  conditions  in  an  iterative  sense.  The  end  plate  is  considered  as  two 
curves,  an  upper  curve  and  a  lower  curve.  The  curves  are  no  longer  defined  as  halves 
of  a  parabola;  instead,  they  are  defined  as  polynomials  which  must  match  points  and 
slopes  at  either  end.  The  expression  for  the  cavity  surface  y(x)  is 

y(x)  =  ax  +  by/ A  -  x  +  c  (C.4) 

where  the  coefficients  a,  b,  c  are  found  from  the  conditions 

y(o)  =  VT 

j/(0)  =  y'j 

y{*)  =  yi • 

yr  and  y'T  are  the  y-position  and  slope  of  the  cavity  at  the  juncture  with  the  end 
plate,  as  shown  in  Figure  C-9.  yi  is  the  y-position  at  the  cavity  trailing  edge.  Notice 
that  equation  C.4  forces  y(x)  to  have  infinite  slope  at  x  =  A. 

At  this  point,  it  was  also  recognized  that  the  extent  of  the  upper  and  lower  curves 
need  not  be  the  same.  There  are  therefore  two  defining  parameters,  Aj  and  Au  (see 
Figure  C-8).  At  the  end  of  every  iteration,  a  new  cavity  shape  is  found  by  integrating 
the  source  strength  forward  from  the  leading  edge  on  the  suction  side  and  from  the 
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Figure  C-7:  Results  of  linear  theory  and  nonlinear  theory  using  the  end  parabola 
model.  N  AC  A 16006  foil  at  a  =  10°  and  a  =  18°.  The  cavity  length  is  £  =  2.0;  the 
cavity  detaches  at  the  leading  edge.  The  length  of  the  transition  zone  is  A  =  .05 


Figure  C-8:  A  sketch  of  the  curved  end-plate  model  for  supercavitation. 

trailing  edge  on  the  pressure  side.  The  shift  of  the  cavity  surface  is  taken  along  its 
normal  everywhere  except  in  the  transition  zone,  where  we  allow  only  the  vertical 
component  of  the  shift.  The  cavity  closure  condition  requires  that  the  vertical  shifts 
A  of  the  upper  and  lower  halves  of  the  end  plate  be  equal  (as  shown  in  Figure  C-9). 
In  this  way,  a  new  cavity  shape  has  been  found  and  is  ready  for  the  next  iteration. 
However,  at  this  point  we  regenerate  the  end  plate  using  equation  C.4.  In  other  words, 
we  force  the  surface  to  be  continuous  with  the  current  iteration.  This  is  repeated  until 
the  nonlinear  solution  converges. 

The  curved  end  plate  model  is  the  most  robust  model  I  tried  (i.e.,  it  allowed 
smooth  convergence  at  very  high  angles  of  attack).  In  a  sense,  this  model  is  a  more 
appropriate  “Riabouchinsky-analogue”  than  the  flat  end  plate  model,  because  it  al¬ 
lows  the  upper  half  of  the  plate  to  be  longer  than  the  lower  half.  It  turns  out  that 
this  is  necessary  in  order  to  obtain  smooth  (continuous  slope)  pressure  distributions 
(corresponding  to  continuous  cavity  curvature)  at  the  cavity /end- plate  juncture.  This 
is  shown  in  Figure  C-10,  where  computations  are  shown  for  a  flat  plate  hydrofoil  at 
a  =  30°  with  two  choices  of  A j  and  Au.  In  the  first  computation,  (Aj,Au)  =  (.22,  .33) 
and  the  resulting  pressures  are  shown  to  overshoot  the  cavity  pressure  in  the  transi¬ 
tion  zone.  In  the  second  computation,  (Aj,  Au)  =  (.13,  .23)  and  the  resulting  pressures 
are  slope-continuous  at  the  juncture  and  everywhere  greater  than  the  cavity  pressure. 
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Figure  C-9:  A  blow-up  of  the  transition  zone,  showing  how  the  cavity  is  updated 
between  iterations  with  the  curved  end  plate  termination  model. 

C.2.2  Pressure  Law 

A  pressure  model,  similar  to  the  one  used  for  partial  cavitation,  was  also  imple¬ 
mented  for  supercavitation.  This  model  is  characterized  by  three  parameters,  Ai ,  A2 
and  as  shown  in  Figure  C-ll.  Ai  and  A2  are  two  adjacent  zones  (unlike  the  end 
plate  model,  they  do  not  differ  from  top  to  bottom).  One  zone  extends  upstream 
a  distance  Aj  from  the  original  cavity  “trailing  edge”  (corresponding  to  the  cavity 
length  assumed  for  the  linear  solution).  The  other  zone  extends  downstream  a  dis¬ 
tance  A2.  The  cavity  is  assumed  to  be  open  and  the  openness  is  specified  to  be  S*. 
The  cavity  closure  condition  is  equivalent  to  the  one  used  for  the  curved  end  plate 
model,  requiring  that  the  upper  and  lower  cavity  coordinates  at  the  upstream  end  of 
the  transition  zone  have  the  same  vertical  translation  from  one  iteration  to  the  next. 
The  wake  is  assumed  to  be  open  and  the  openness  extends  to  downstream  infinity. 

The  effect  of  the  pressure  law  parameters  on  the  cavitation  number,  cavity  volume 
and  lift  and  drag  coefficients  is  shown  in  Table  C.l.  Notice  that  the  length  of  the 
transition  zones  A2,A2  have  a  fairly  large  effect  on  <r  and  V,  while  the  forces  are 
relatively  independent  of  the  parameters. 

A  direct  comparison  between  the  curved  end  plate  model  and  the  pressure  law 
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Figure  C-10:  Cavity  shapes  (from  linear  theory  and  nonlinear  theory)  and  pressure 
distributions  (from  nonlinear  theory)  for  a  flat  plate  hydrofoil  ata  =  30°.  The  slope 
continuity  of  the  pressure  distribution  is  shown  to  depend  on  the  parameters  Xi  and 


Figure  C-ll:  A  depiction  of  the  pressure  law  for  supercavitation. 
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Table  C.l:  The  effect  of  parameters  on  the  cavity  solution  using  the  pressure  law  ter¬ 
mination  model.  Computations  are  for  a  N  AC  A 16005  thickness  form  with  a  NACAa.8 
camber  profile  with  a  maximum  camber  of  .03.  a  =  4°,  j  =  1.5,  ^  =  0.01  and  A  =  0.1. 


Table  C.2:  The  difference  between  the  curved  end  plate  model  and  the  pressure 
law  model  for  the  same  flow  geometry  as  in  the  previous  table.  The  pressure  law 
termination  model  parameters  are:  \\  —  0.1,  A2  =  0.3,  6*  =  0.1,  and  A  =  0.1.  The 
end  plate  termination  model  parameters  are:  Aj  =  0.11  and  Au  =  0.15. 

model  is  shown  in  Table  C.2  for  the  same  flow  configuration  treated  in  Table  C.l. 
Shown  also  with  the  table  are  the  computed  cavity  shapes,  showing  the  difference 
termination  models. 
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